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ABSTRACT 

m ' 

1 We perform 3-dimcnsional hydrodynamic simulations of gas flowing around a plan- 

etary core of mass M p i=l 0M ffi embedded in a near Keplerian background flow, using 
a modified shearing box approximation. We assume an ideal gas behavior following 
an equation of state with a fixed ratio of the specific heats, 7 = f .42, consistent with 
the conditions of a moderate temperature background disk with solar composition. 
■ No radiative heating or cooling is included in the models. We employ a nested grid 

hydrodynamic code implementing the 'Piecewise Parabolic Method' with as many as 
six fixed nested grids, providing spatial resolution on the finest grid comparable to the 
£^ ■ present day diameters of Neptune and Uranus. 

I We find that a strongly dynamically active flow develops such that no static en- 

CS . velope can form. The activity is not sensitive to plausible variations in the rotation 

curve of the underlying disk. It is sensitive to the thermodynamic treatment of the 
gas, as modeled by prescribed equations of state (either 'locally isothermal' or 'locally 
isentropic') and the temperature of the background disk material. The activity is also 
' sensitive to the shape and depth of the core's gravitational potential, through its mass 

, and gravitational softening coefficient. Each of these factors influence the magnitude 

■<^j- ■ and character of hydrodynamic feedback of the small scale flow on the background, 

' and we conclude that accurate modeling of such feedback is critical to a complete 

understanding of the core accretion process. 

The varying flow pattern gives rise to large, irregular eruptions of matter from 
I the region around the core which return matter to the background flow: mass in the 

envelope at one time may not be found in the envelope at any later time. No net mass 
accretion into the envelope is observed over the course of the simulation and none is 
expected, due to our neglect of cooling. Except in cases of very rapid cooling however, 
as defined by locally isothermal or isentropic treatments, any cooling that does affect 
j_h " the envelope material will have limited consequences for the dynamics, since the flow 

quickly carries cooled material out of the core's environment entirely. The angular 
momentum of material in the envelope, relative to the core, varies both in magnitude 
and in sign on time scales of days to months near the core and on time scales a few 
years at distances comparable to the Hill radius. The dynamical activity contrasts 
with the largely static behavior typically assumed within the framework of the core 
accretion model for Jovian planet formation. 

We show that material entering the dynamically active environment may suffer 
intense heating and cooling events the durations of which are as short as a few hours 
to a few days. Shorter durations are not observable in our work due to the limits 
of our resolution. Peak temperatures in these events range from T ~ 1000 K to as 
high as T ~ 3 — 4000 K, with densities p ~ I0~ 9 — I0~ 8 g/cm 3 . These time scales, 
densities and temperatures span a range consistent with those required for chondrulc 
formation in the nebular shock model. We therefore propose that dynamical activity 
in the Jovian planet formation environment could be responsible for the production 
of chondrules and other annealed silicates in the solar nebula. 
Key words: hydrodynamics, planets: formation, planetary systems: formation 



© 0000 RAS 



2 Andrew F. Nelson and Maximilian Ruffert 



1 INTRODUCTION 

The classical 'core accretio n' model for Jovian planet forma- 
tion IjWuchterl et fflZj|200CO is characterized by the idea that 
a 5-15 earth mass (M©) rock/ice core forms at a distance 
of ~5 AU or more from the central star (|Bosdl 19951 ). with a 
gaseous envelope surrounding it. At some critical point dur- 
ing the growth, the gaseous envelope becomes unstable and 
collapses, ultimately accreting as much as several times the 
mass of Jupiter. 

In outlin e, the core accretion model consists of three dis- 
tinct phases (|Pollack et al ] |l996h . First, a sold core accretes 
from smaller solid bodies like comets and asteroids in the 
solar nebula. The mass of this core may be as much as 20-30 
M®, depending on the model, though values of 5-15M® are 
much more typically expected. The relatively rapid growth 
of this core is cut off when it depletes its 'feeding zone', at 
which point begins a much slower period of growth of both 
the core, from material migrating into the feeding zone from 
other parts of the disk, and of the surrounding hydrostatic 
gaseous envelope. In this second phase, the gas accretion 
rate is regulated by the balance of radiative cooling and the 
continued heating of the envelope by a low rate of solid body 
accretion. When the combined mass of the core and envelope 
reaches a (model dependent) critical value of ~ 30M®, the 
hydrostatic envelope becomes unstable and period of rapid 
gas accretion follows, defining the third, 'core instability', 
stage of evolution. The planet's final mass is determined 
during this phase. 

The most critical problem afflicting the core accretion 
model is that the sec ond stage may requ ire as much as 6 
million years or more (|Pollack et alj|l996h . but the circum- 
stellar disk out of which the planet pur portedly forms ma y 
have a lifetime of only 2-4 million years l|Haisch et aZ.|[200ll ) . 
More critically, when the dynamical consequences of the in- 
teraction between the forming planet and disk are taken 
into account, the planet's lifetime may shrink to only a few 
x 10 4 years before it mig rates inwards toward the star and is 
accreted by it (see e.g. IWard &: Hahnl l2000. and references 
therein). While some migration is certainly necessary, since 
the vast majority of extra solar planets so far discovered have 
been observed in orbits less than (often much less than) the 
^5 AU where they are expected to form, the rates expected 
from theory would imply that no planets could outlive the 
disks in which they form. 

1.1 Newer wrinkles in core accretion 

IPollack et all jl996) show that among the most sensitive 
input parameters for the formation time scale is how the 
planctesimal accretion rate varies with time, especially in 
the second phase of growth after the core has grown sub- 
stantially and depleted its feeding zone. The cause of the 
sensitivity is the continued heating of the envelope gas by 
the solids as they pass through and are dissolved in the enve- 
lope when they accrete. The heating then stalls both further 
envelope contraction and gas accretion onto the core. 

These models were performed in the context of a sin- 
gle, isolated core in a background disk where no other cores 
competed for the same resources and no migration through 
the disk occurred. The time scale can be s hortened con- 
siderably when either assumption is relaxed. lAlibert et all 



have shown that including migration can shorten the 
required lifetime to < lMyr even in a much less massive 
disk than originally studied because, due to the migration 
the planet does not deplete its feeding zone for solid bodies. 
Instead, its heavy element core continues to grow rapidly, 
triggering runaway gas accretion much more quickly than 
otherwise. A possibly undesirable consequence of the model 
is that critical core mass is higher, perhaps in conflict with 

the data for the observed core mass of Jupiter. 

Although not including migration, Irlubickvi et all 

(2005) have simulated the effects of competing embryos with 
a 'cutoff' of the solid body accretion when the core reaches 
a predetermined mass, modeling the fact that the solids at 
more than some critical distance from one core have been 
accreted by another. In addition, they study the sensitivity 
of their models to the treatment of opacity in the envelope, 
which governs the rate at which cooling of the envelope (fol- 
lowed by contraction and further accretion) can proceed. 
Previous studies have used opacities that are appropriate 
for interstellar grains and therefore do not properly model 
the processing expected during the evolution of the circum- 
stellar disk and later the gaseous planetary envelope itself. 
Such processing, especially in the envelope, means that opac- 
ities will be reduced as grains coag ulate and 'rain out' o f 
the forming planet's upper envelope. iHubickvi efai\ (I2005T ) 
model such processing with a temperature dependent reduc- 
tion of the interstellar opacity value and, coupled with the 
cutoff mass, find that the required formation time scale is 
reduced to as little as 1 Myr, with critical core masses of 
5-10M®. 

1.2 Shortcomings in the models and goals for this 
work 

An important underlying assumption made in core accre- 
tion models is that a hydrostatic envelope structure exists, 
so that a one dimensional approximation may be used to re- 
produce the evolution. On the other hand, previous numeri- 
cal studies of planets embedded in disks have shown that this 
hydrostatic assumption may n ot be realistic, even fo r very 
low mass cores. For example, iD'Angelo et all (|2002) show 
that in 2D nested grid simulations, spiral structure exists 
in the circumplanetary disks of cores as small as 3M®, and 
that for cores of > lOM® strong spiral shocks are present. 
iBate et~al\ (|2003l ) perform similar simulations in 3D and also 
observe spiral structures and rotation, though they do not 
find shock structures as strong as in the 2D work. Their 
simulations implemented an isothermal equation of state 
however, which is certainly inappropriate for such dense re- 
gions as the envelope of a forming planet and may result 
in flows that do not reproduce the t r ue pic ture of the evo- 
lution. Later work of lAvliffe fc Bate! (|2009l . hereafer AB09) 
relaxes the isothermal assumption in some of their simula- 
tions, instead using a flux limited radiative transfer model of 
a small cutout region of the disk whose evolution was simu- 
lated with a Smoothed Particle Hydrodynamics (SPH) code 
with a resolution of 3 x 10 4 SPH particles. Their conclu- 
sions are similar to the isothermal results in the sense that 
in both cases, spiral structures and/or envelope asymmetries 
develop, but are weaker in the radiative transfer models, and 
accretion rates onto the core are lower. 

Although these studies are good first steps, to date 
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there has been little exploration of the region of parame- 
ter space between the two extremes of, at one end, detailed 
modeling of the disk with very coarse modeling of the core 
and envelope, and at the other, of detailed modeling of the 
core and envelope with coarse or non-existent modeling of 
the disk. The simulations of AB09 for example, which are 
among those with the highest claimed resolution of the re- 
gion around the core, employ a total of only 3 x 10 4 or 
1 x 10 s (in different simulations) particles over the entire 
volume. In order to fully understand the core accretion sce- 
nario, detailed models of the flow around and onto the core 
at much higher resolution are necessary. 

Models of accre tion ont o com p act o bjects have been 
published (see e.g. iRuffertJ 1 19971 . Il999l . and references 
therein), an d in many cases large sca le instabilities in the 
flow develop iFoglizzo fc Ruffertl 1 19991 ) that are often char- 
acterized by 'flip-flop' behavior in the sense that angular 
momentum accretion can temporarily switch sign. While 3D 
models produce uniformly less active flow than 2D models, 
both produce such behavior and in most cases, the accretion 
rate of various flow variables is comparable to the classical 
Bondi-Hoyle rate even in flows with transverse gradients. 

The applicability of these models and their results to 
planet formation would appear to be quite limited however 
because they typically assume a central body which accretes 
matter with perfect efficiency (e.g. a 'black hole'). In the con- 
text of planet formation, the accretor is of course a rocky 
planetary core, onto which gas accretion is forbidden and 
an envelope structure develops. Moreover, while the back- 
ground gas flow in previous studies may include some trans- 
verse velocity or density gradients, it is uni-directional. In 
contrast, the flow around a forming planet is quite complex, 
with a flow that circulates in opposite directions relative to 
the planet inside and outside its orbit, but which librates in 
a region near its orbit. The existence of the planet may drive 
complex interactions between each of the three regions and 
with planet's envelope as well. 

With such large differences in the accretion assumptions 
and in the background flow pattern, will similar sorts of 
dynamical activity as in the simpler cases occur? How does 
the existence of dynamical activity (if indeed it is present) 
affect the time scale required for the core accretion model 
to succeed? In other words, can dynamical activity somehow 
short circuit the slow, near hydrostatic evolutionary process 
assumed in the second phase of planet growth, replacing it 
with a more rapid growth phase that avoids the problems 
of rapid Type I migration? Or does it instead lengthen the 
process? 

The goal of this paper is to begin an investigation of 
the accretion process in order to determine the character of 
any dynamical activity that may occur. We will extend the 
previous comparatively low resolution, two and three dimen- 
sional calculations exploring planet migration to very high 
resolution three dimensional calculations specifically study- 
ing the mass accreting onto a IOMq core. The models in our 
study will consist of a small cubic 'cutout' region of a cir- 
cumstellar disk and centered on the core. They will therefore 
not include migration of the planet itself, but will include 
some effects of the background on the accretion flow, such 
as the effects of different disk density and temperature dis- 
tributions. In section [2] we define the physical conditions 
assumed in our models and numerical methods used to real- 



ize the evolution, as well as the initial conditions we imple- 
ment for our numerical models. Section [3] contains the basic 
results of our simulations and introduces some metrics by 
which those results may be compared more quantitatively. 
Section [4] expands on these analyses with further studies of 
the sensitivity of the results to various parameters. Finally, 
section [5] proposes a potential link between the dynamical 
activity around a forming Jovian planet and the formation 
of chondrules. Our concluding remarks are found in section 
[6] where we summarize our conclusions and the implications 
they have for Jovian planet formation and compare our re- 
sults with previous work. Finally, in section [JJ we discuss 
a number of questions that may be profitably addressed by 
future work similar to our own. 



2 THE PHYSICAL MODEL 

In this section we describe the physical models we use to 
study the flow pattern of gas around the forming planetary 
core, the numerical code used to evolve the models forward 
in time and the initial conditions of those models. Our study 
will include or exclude a variety of physical properties of the 
system in different simulations. Therefore we will outline the 
most physically inclusive models, and in later discussion note 
specifically which properties have been included in which 
simulation. In general, the properties we study are either 
substitutive, eg. either one equation of state or another, or 
cumulative, e.g. gravity from different components of the 
system, so that separate, excluded components can simply 
be neglected in the specification of the initial state. 



2.1 The mathematical framework, numerical 
realization and physical assumptions 

In previous works studying disks, it has been customary to 
implement either a fully curvilinear coordinate system (usu- 
ally cylindrical) or a 'shearing sheet' approximation. The 
latter form is an approximation in which locally Cartesian 
coordinate system is used rather than the full equations for 
curvilinear coordinates which simplifies the analysis, espe- 
cially in the context of linearized instability analyses. It ne- 
glects most of the non-uniform geometric features of the 
coordinate system and instead adds a number of terms to 
model approximately some of its effects, and the rotation of 
the coordinate frame itself. 

We have implemented neither procedure here, but 
rather an intermediate step between them which we call a 
'modified shearing sheet' coordinate system. As does the 
standard shearing sheet coordinate system, this modified 
shearing sheet system neglects the geometric terms in the 
full curvilinear system of equations. Its main difference from 
the standard form is that it retains a more precise treatment 
of the gravity and fictitious forces on the gas, as noted below. 
Due to our somewhat unusual treatment of the coordinate 
system, we shall describe here in some detail both the deriva- 
tion of the background analytic and numerical frameworks, 
as well as the initial and boundary conditions used in our 
models. 
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2.1.1 Mapping cylindrical to modified shearing sheet 
coordinates 

Like the shearing sheet, our mathematical formalism re- 
quires that the true cylindrical coordinate system be 
mapped onto a Cartesian coordinate system centered on the 
planet. We define the mapping so that the three pairs (x, 
r), (y, r<p) and (z, z) each correspond as follows. Mapping 
the cylindrical coordinates of the disk to the Cartesian co- 
ordinates of our grid requires the three identifications 



x — r — dpi; y = r(f>; z = z 

with dpi defined as the semi-major axis of the planet. This 
identification means that each x grid coordinate is defined 
to be at the same distance from the star even though this 
condition would not be true for a true Cartesian grid cen- 
tered on the planet. The errors that are made by such an 
identification will be small as long as the extent of the co- 
ordinate grid itself is small, compared to the distance to the 
origin. 

In our models, we typically simulate the gas flow in 
a region ~ 0.45 AU (4 'Hill radii' from a 10M® object) 
in each coordinate direction from a point centered 5.2 AU 
away from the origin. In a cylindrical coordinate system, 
two 'cubic' volume elements at the inner and outer edges 
of a grid of this size (i.e. rdrd(j>dz, with each side of the 
volume identical in length) will differ in actual volume by 
a factor ~ 1.4, whereas in our Cartesian coordinate system 
they will be equal. Because the flow is largely azimuthally 
directed, very little gas will undergo (or fail to undergo in 
our coordinate system) expansion or contraction by such a 
factor and therefore only small errors will develop from this 



2.1.2 The equations of hydrodynamics 

Here we state explicitly the fully expanded form of the equa- 
tions we solve. The hydrodynamic variables for mass, mo- 
mentum and energy conservation are defined with the usual 
symbols with p referring to the mass density, v x , v y and v z 
to the velocities in each of the three coordinate directions 
and E = p{u + v 2 /2) to the total energy, where v is the 
magnitude of the total velocity and u is the specific internal 
energy. 

The continuity equation for mass conservation is given 
by 



dp d(pv x ) d(pvy) d(pv z ) = Q 
dt dx dy dz 



(2) 



and the equations of motion in the three Cartesian directions 
are given by 



d(pv x ) d(pv x v x ) d(pv x v y ) d(pv x v z ) 



PVy 2 



dt 



dx 



9y 



dz a p i + x 

dp d® 



dx 



dx 



(3) 



d(pvy) d(pv y v x ) d(pVyVy) d(pv y v z ) pv x v y 



dt 



+ 



dx 



+ 



dy 



+ ■ 



dz 



+ 



a p i + x 



+2 



PV X Vy 



dp d$ 

= — pit ( 4 ) 

a p i +x dy dy 



and 

d{pv z ) ^ d(pv z v x ) d(pv z v y ) d(pvzVz) 



dt 



dx 



dy 



dz 



dp d$ 
= -d~z- p ^- (5) 
For models in which an equation for the conservation of en- 
ergy is required (i.e. models without locally isothermal or 
locally isentropic equations of state), the energy conserva- 
tion equation is 



(1) dE | d(v x (E + p)) | d(v v (E + p)) | d(v z (E + p)) 

dt dx dy dz 

In equations [3] and U above, V y is defined as 

Vy = Vy + (dpi + X)tl p l (7) 

where fi p i is the orbital angular velocity of the planet. V y 
corresponds to the total azimuth velocity in a non-rotating 
cylindrical coordinate frame. An equation of state for the gas 
closes the system of equations and is described in section 
12.1.31 The form of the total gravitational potential, $, is 
discussed in section [2.1.41 below. 

Note that the above set of equations include no explicit 
contribution due to any viscous properties of the fluid. Such 
effects are well known to be present in the background cir- 
cumstellar disk, with their origin commonly ascribed to tur- 
bulence g enerated by the so-called Magn eto-Rotational in- 
stability (jBalbus. Hawlev fc Stondll996l ). Both the spatial 
and time scales relevant for this source of dissipation are 
inapplicable to the immediate environment of a core em- 
bedded in a circumstellar disk however, and we therefore 
neglect turbulence as a source of dissipation in our simu- 
lations. We similarly neglect other, less specific, sources of 
dissipation in favor of modeling only those characteristics of 
the flow to which we can define explicitly. Coupled with our 
numerical code (section [2TT3}, which itself contributes little 
dissipation to its solution of the flow equations, our simula- 
tions will therefore be highly inviscid. The only significant 
sources of dissipation will originate in shocks that develop 
in the flow, which the numerical scheme is known to handle 
well. 

For reference, we point out the differences between our 
treatment and both the fully curvi-linear coordinate system 
and the shearing sheet, in Appendix [X] 



2.1.3 The thermodynamic treatment of the gas 

The system of hydrodynamic equations are closed assuming 
an ideal gas equation of state given by: 



(7 - l)pu 



(8) 



where 7 is the ratio of specific heats and p, p and u are re- 
spectively the pressure, density and specific internal energy 
of the gas. In order to derive the global initial condition 
for our simulations, we employ a vertically integrated (2D) 
equation of state of the form: 



P = (7 - l)£u 



(9) 



where P and E are the vertically integrated pressure and 
the mass surface density respectively. This form is exactly 
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analogous to that used in 3D for the simulations themselves, 
so that no adjustments to the rotation curve or any other 
quantities are required to establish a steady state initial con- 
dition. Finally, the sound speed and temperature are related 
by the equation 



(10) 



where R is the gas constant and \i = 2.31 is the mean molec- 
ular weight of a gas of solar composition. 

The choice of 7 is not obvious in the range of temper- 
atures and densities we study because in various regimes, 
some or all of a number of internal molecular states may 
be important. In order to avoid such complications, we as- 
sume that the gas is of solar composition with an effective 
single component accounting for both hydrogen and helium, 
so that the value for 7 = 1.42. This value implies that the 
rotational degrees of freedom in hydrogen gas are active (cor- 
responding to gas temperatures > 100 K) but its vibrational 
degrees of freedom are inactive (corresponding to tempera- 
tures < 800 — 1000 K), and will be most representative of 
moderate temperature regions of the disk expected through- 
out most regions of our simulations. 

Apart from heating from hydrodynamic processes like 
PdV work and shocks, we impose no other assumptions re- 
garding the heating and cooling processes that are active 
in the disk. Of particular importance in this regard is our 
neglect of radiative heating and cooling, which will undoubt- 
edly modify the behavior of the models to some extent. As 
a first, crude attempt to model such effects, we also perform 
a number of simulations using either a locally isothermal or 
locally isentropic equation of state, defined respectively as 



P = pc s 
and 

p = Kp< 



(11) 



(12) 



with the value of 7 set to 1.01 or 1.42 for the locally isother- 
mal or locally isentropic cases respectively. The physical in- 
terpretation these equations of state in terms of their impli- 
cation s for heating and co o ling a re discussed in I Nelson et all 
l|2000h and iPickett et all (|2003l '). The values of c\ and K 
are set by the initial conditions of the simulation and are 
thereafter fixed for its duration. As with the ideal gas case, 
vertically integrated forms are employed to determine the 
global initial condition, prior to the beginning of the actual 
simulation. 



2.1.4 The treatment of gravity 

Gravity from the central star, from the planet core and from 
the circumstellar disk are included in our simulations. The 
contribution from the disk is further split into two parts, 
representing the portion of the disk inside and outside the 
simulation box. 

The total gravitational potential is defined as a sum of 
several components from the star $», the planet <3> p i, the 
background disk $0, and the mass inside our simulation 
cube $i oc as: 



The last term, $i oc , describing the initial value of the lo- 
cal gravitational potential due to material in the simulation 
cube, is required to avoid double counting disk mass that is 
otherwise computed in both the background disk and then 
again locally in the simulation cube, as we discuss in more 
detail below. 

We assume that the central star is fixed at all times 
during our simulations, which means that its contribution 
to the gravitational force on the gas is constant. The gravity 
from the planetary core is also fixed in time, relative to the 
grid. Both the star and the planet are modeled as softened 
point masses with a softeni ng that is similar in form to a 
Plummer law (|Ruffertll 19941 ) with a potential given by 



$ = 



—GM 



(r 2 + e 2 <5 2 exp[-r 2 /e 2 5 2 ]) 1 / 2 



(14) 



$ = $, + $ P 1 + $D + $ 



(13) 



where G is the gravitational constant and the quantities $ 
and M are replaced by with variants subscripted for the star 
or the planet planet as appropriate. The softening parameter 
e is a tunable constant multiplying the grid spacing 8, and 
in our simulations is set to unity. The exponential term in 
the denominator causes the softening to decrease to zero at 
distances of more than a few grid spacings from the point 
mass, where it is not needed. 

The distances, r, between the star or planet and each 
individual grid cell are determined from the grid spacing it- 
self, with one important modification. In order to simplify 
our initial condition (see section [2.2.2|) . we neglect the y and 
z components of stellar gravity, so that the distance from 
the star to any point in the grid includes only the x compo- 
nent, effectively treating the x coordinate (plus the offset, 
dpi, defining the planet's semi-major axis in the underly- 
ing circumstellar disk) as a radial coordinate. The planet's 
contribution includes all three Cartesian components of the 
force. 

The gravitational potential of the disk is split into two 
components, one from the disk as a whole (the 'global' 
component), and a second from the portion of the disk in- 
cluded in our cubic cutout region (the 'local' component). 
The global disk self gravity is computed once for the ini- 
tial state, and thereafter remains constant for the duration 
of the simulation. The local potential is dependent on the 
time dependent flow through the disk and must therefore be 
recomputed at each time step. 

We calculate the global component from the specified 
global mass distribution of the disk (see 12.2.11 below) using 
a Fourier transform based gravi tational potential solver i n 
cylindrical coordinates (see e.g. iBinnev fc Tremaineill987l ). 
Although the disk is axisymmetric, so that we would expect 
that a one dimensional solution for the potential with such 
a solver would be adequate, we have found that in practice 
this solver requires that a number of azimuth cells be used 
in order to converge to a well determined potential value. As 
discussed in sect ion 12. 2 . 1 1 b elow . the radial grid for our global 
initial condition requires as many as approximately 30000 
cells for our highest resolution models. We have found that 
with this radial resolution we require at least 1024 azimuth 
grid cells be included in the potential solver to produce a 
potential solution that is accurate to numerical double pre- 
cision (i.e. one part in ~ 10 15 ) over the entire radial range 
of 0.5 to 20 AU included in the global initial condition. Al- 
though we have made no specific investigation of the reasons 
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for the required number of azimuth cells, it seems likely to 
be due to the existence of a small self-contribution to the 
potential as obtained from the Fourier t ransform technique 
in cylindrical coordinates (|Massetlll997h , for which we at- 
tempted no correction. As for the stellar gravity, the back- 
ground disk gravity is mapped to the x coordinate only. 

The gravitational potential of material on the grid, $i oc , 
is determined from the Poisson equation 



8 



d 



d 



dx 2 dy 



dz 



(15) 



We compute its value at each time using a 3D Cartesian co- 
ordinate FFT based potential solver. Due to the nested grid 
structure and the inclusion of the global disk potential there 
are two additional complications regarding the calculation 
of the potential that do not occur in the standard method 
of solution on a single grid. Both complications can be cir- 
cumvented by straightforward application of the fact that 
the Poisson equation is linear, so that separate terms may 
be added and subtracted independently. 

The first complication is that the potential calculation 
on each of the nested grids require contributions from mass 
outside of the sub-cub e on which the potential is to be cal- 
culated. iRuffertl l|l992tl , notes that the full potential can be 
obtained from a sum of the contribution from a given fine 
grid, and from a second potential calculation that is per- 
formed on the overlying coarse grid, in which the mass on 
the finer grid is temporarily deleted. The portion of this sec- 
ond potential (on a given coarse grid) that overlaps the fine 
grid is mapped onto that fine grid as a background term and 
summed with the fine grid computation to obtain a correct 
final value. 

The second complication arises because the global grav- 
itational potential of the disk accounts for the contribution 
due to mass located both inside and outside of our local sim- 
ulation box. This means that simply adding the contribution 
from the mass distribution within the simulation box would 
result in that matter being accounted for twice in the net 
gravitational potential. To determine the correct local po- 
tential, we first compute the gravitational potential of mat- 
ter in the box due to the mass distribution of the initial state 
and save this value for the life of the simulation. The initial 
contribution is subtracted from the potential calculated at 
all later times, so that effectively only the difference in the 
mass distribution of the later and initial states is accounted 
for in the local potential and the double counting is avoided. 

2.1.5 The numerical code 

To model the flow at the extremely hi gh spatial reso lution 
required we have modified the code of iRufferj (|1992r h This 
code is based on the 'Piecew ise Parabolic Method' (PPM) of 
IColella fc Woodward! (|l984l ). in which a high order polyno- 
mial interpolation is used to determine cell edge values used 
in calculating a second order solution to a Riemann prob- 
lem at each cell boundary. The interpolation is modified in 
regions of sharp discontinuities to track shocks and contact 
discontinuities more closely and retain their sharpness, while 
a monotonizing condition smoothes out unphysical oscilla- 
tions. The solution to the one-dimensional Riemann problem 
is then used to calculate fluxes and advance the solution in 
time. Since the Riemann problem solution explicitly models 



the physical dissipation present in shock structures, no arti- 
ficial viscosity is required for stability of the code, and none 
is included. 

This code solves the hydrodynamic equations on a three 
dimensional Cartesian grid, and includes the possibility of 
including a series of statically generated nested grids. Each 
grid contains an identical number of zones in each of the 
three coordinate directions, but the grid spacing decreases 
by a factor of two when proceeding from a coarse to fine 
grid. 

We have modified the original, strictly Cartesian code 
to account for centrifugal and Coriolis forces required in a 
coordinate system rotating with the planet core, but have 
neglected terms required to reproduce a fully curvi-linear 
coordinate system, as noted in section [2.1.2l above. This ap- 
proximation will remain valid so long as the region modeled 
is relatively small, so that the curvature is unimportant. 

Depending on the simulation (see table [T] below) , be- 
tween three and six nested grids each with 64 3 zones are 
used to model a cubic region around the planet. Successively 
finer grids have the same number of zones, but half of the lin- 
ear dimensions of the immediately coarser grid and are also 
centered on the planet. When the simulation volume is de- 
fined by a region ±4 Hill radii (section l2,1.6[) around a lOM® 
planet, the linear resolution on the finest grid of our high- 
est resolution models (6 nested grids) is 8x ~ 6.5 x 10 9 cm 
in size, about 30% larger than the diameter of Neptune or 
Uranus and the spatial volume is resolved with a total of 
approximately 1.4 million zones. 



2.1.6 The Hill and accretion radii 

In the discussion throughout this paper, two quantities will 
appear repeatedly, and so we introduce them here. They are 
the accretion radius of the planet core, defined as 



Ra 



GM„ 



(16) 



and its Hill radius 



(17) 



where G is the gravitational constant, M v \ and M* are the 
mass of the planet core and the central star, respectively, 
dpi is the semi-major axis at which the core orbits the star 
and c 2 is the square of the sound speed of the gas in the 
unperturbed (initial) disk near the core. 

Taken in its original context of the Bondi accretion 
model, the accretion radius defines the distance from the 
core at which the magnitude of the internal energy of an 
isothermal gas and its gravitational potential energy become 
comparable. It therefore represents a measure of the spatial 
scale over which exchanges between the fluid's thermal en- 
ergy and energy bound up in the bulk kinetic motion derived 
from its infall onto the core might be expected to lead to 
dynamical activity. Together with the mass of the core, the 
background temperature, through its presence in the rela- 
tion for the sound speed, completely specifies the accretion 
radius. 

The Hill radius defines another spatial scale of inter- 
est. Its value is derived from the Jacobi integral, a quantity 
which plays a role similar to the total energy in other (in- 
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ertial) systems, and which is the only conserved quantity 
in the restricted three body problem. It relates the gravita- 
tional potential energies imposed by the core and star on a 
test particle with that particle's kinetic energy as it travels 
through the combined system and a fourth 'centrifugal po- 
tential' term quantifying the effects of the non-inertial refer- 
ence frame assumed for the system. It also defines the inner 
and outer extents of the Roche lobe of the planet, which is 
the full 3D volume around the core for which the combined 
gravitational influences of the core and star are compara- 
ble to the influence of the fictitious forces arising from the 
rotating reference frame. 

2.2 Initial conditions 

We develop the initial conditions for our simulations in a 
two stage process. First, we define a set of global conditions 
that describe an entire circumstellar disk in two spatial di- 
mensions (r,(j>). From these conditions, we extract a three 
dimensional cubic 'cutout' region centered on a point mass, 
assumed to be embedded within the disk. This dual speci- 
fication requires a small note of clarification for readers, in 
that in various places throughout the text, 'mass density' 
might refer to either a mass per unit area, or to a mass per 
unit volume. While in most contexts the distinction will be 
clear, we will typically refer to each as either 'surface' or 
'volume' density, respectively. 

2.2.1 The Global Model 

The initial conditions are developed from a global model of 
an accretio n disk quite sim i lar in morphology to the disks 
modeled in INelson fe Bend {2003). Specifically, we assume 
that a point mass representing the central star is surrounded 
by a circumstellar disk, whose inner and outer radial dimen- 
sions are E4 = 0.5 AU and R a = 20 AU. 

The mass in the disk is distributed according to a sur- 
face density power law of the form 

EM = £ P i(^) P , (18) 

with a similar power law for temperature taking the form 

T(r) = T pl (^i) ? . (19) 

E p i and T p \ are the surface density and temperature at the 
assumed orbit radius, a p \ of the planet, respectively. The 
power law exponents are ordinarily set to p — q — 1.0. In 
order to study the sensitivity of our results to the mass and 
temperature distributions however, we have also performed 
simulations with other exponents as well (see table HJ. 

The surface density coefficient is set to E p i = 
500 g/ cm 2 . This value is som ewhat smaller than that implied 
by the lPollack et all (|l996h model, for which 6 Myr were re- 
quired to grow a Jovian mass planet. It is chosen so as to be 
consistent with the minimum mass solar nebula (MMSN), 
given the dimensions of the disk we assume. With the given 
spatial dimensions of the disk this surface density implies a 
disk mass of Md=0.0358A/q , comparable to the traditional 
MMSN value, but well below the M D ~ 0.1-0.14Mq requi red 
for the classical core accretion model |Pollack et al 1 ll996l ) to 
produce a fully formed Jovian planet within ~ 6Myr. 



The temperature exponent value of q = 1 is chosen to 
ensure that the scale height of the disk is approximately con- 
stant. It is steeper than the value usu ally expected from ob - 
servations of circumstellar disks (e.g. iBeckwith et al.|[l990l ). 
but over the limited radial range we simulate, we do not 
expect the temperature gradient to play a large role. The 
principle benefit of choosing this value will be to simplify the 
correspondence between the global and local disk models, as 
discussed below. The temperature scale is set to T p \ = 200 K, 
comparable to (but slightly warmer than) the ice condensa- 
tion temperature expected to be important for the rapid co- 
agulation of grains that begins the planet formation process. 
With this background temperature, the derived isothermal 
scale height of the disk, H = c 3 /Q, normalized by the dis- 
tance from the star is H/r — 0.08, as shown in figure [T] 

We map each of the thermodynamic quantities onto a 
radial grid in cylindrical coordinates, equally spaced in ln(r). 
In order to ensure the most precise initial condition for each 
sub-grid, the local initial conditions are derived separately 
for each sub-grid. The grid spacing of the ID global grid is 
adjusted at each derivation to be equal to the spacing of the 
current local grid. This means that over the full extent of 
the disk (0.5 to 20 AU) we require as many as 30000 radial 
grid points to determine the global initial condition. Since 
the global initial condition is determined only at the begin- 
ning of the simulation, used to determine a corresponding 
local condition and thereafter discarded, such high resolu- 
tion poses a negligible cost in terms of the total time required 
for the simulation. 

From the predefined temperature, surface density and 
equation of state, we determine a pressure gradient for the 
disk. The pressure gradient is then used in conjunction with 
the gravitational forces from the star and disk to determine 
a rotational velocity, under the assumption that radial ac- 
celerations due to pressure gradients, stellar and global disk 
self gravity and centrifugal forces exactly cancel each other 
and no radial motion exists. No initial perturbations are as- 
sumed in any of the models. 

Finally, we combine the various quantities from our 
specification of the global disk conditions together to cal- 
culate a value of the Toomre Q parameter, defined as: 



where k is the local epicycli c frequency and c s is the sound 
speed. Theoretically (see e.g. lBinnev fc Tremaine|[l987l ). val- 
ues of Q < 1 define the conditions for which self gravitating, 
axisymmetric instabilities in disks are un stable to linear per- 
turbations. Numer ical simulations (e.g. INelson et aZj|l998l : 
iPickett et all\l99$ ) have shown that self gravitating spiral 
instabilities can still grow for values Q < 2 — 3, but that disks 
with higher values are stable. 

The conditions assumed for most of our simulations en- 
sure the global stability of the disk to large scale self grav- 
itating disturbances, as measured by the value of its Q pa- 
rameter. However, in several simulations we will also study 
disks with background temperatures as low as T p \ =50 K. 
Since, in the context of a global initial condition, the value 
of Q in such a disk would certainly fall to values low enough 
(in the outer disk) to allow instabilities to grow so that other 
planet formation mechanisms to proceed, such values may 
seem rather unphysical or at least irrelevant. The simula- 
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Figure 1. Global initial conditions for the disks underlying our simulations. Clockwise from the top left, the four panels of the image 
are the disk surface density, the temperature, the Toomre Q value and the dimensionless scale height. The heavy portion of each line 
corresponds to the portion of the initial condition considered in deriving the local initial condition. 



tions remain interesting however because of their importance 
in outlining the behavior of the flow as a function of the ra- 
tio between the Hill and accretion radii, which we will find is 
an important measure of the dynamical activity in the flow. 

A graphical summary of the global initial conditions 
for our prototype model (discussed in section [3Tl below) is 
shown in figure [T] As indicated, the lower right panel shows 
the Toomre stability parameter, Q, value for the simula- 
tions at each radius. Because Q > 5 for all radii, we expect 
that although global self gravity is included, self gravitat- 
ing disk instabilities would not develop in the background 
flow, validating our assumption that it remains smooth. Disk 
self gravity remains important to consider however, first be- 
cause, from a global perspective, it shifts the position of the 
corotation resonance relative to the planet, and second, from 
a local perspective, because we do not wish to discount the 
possibility that self gravitating instabilities could develop in 
the forming planet's envelope, even if they would not in the 
disk as a whole. 



2.2.2 The Local Model 

We use the global initial condition to derive the local ini- 
tial condition inside our simulation cube in the following 
manner. First, we assume a planet core (modeled as a point 
mass) is embedded in the disk and orbits the star at a dis- 
tance of dpi = 5.2 AU. Its orbit velocity is determined from 
the assumption that its orbit is circular and that it is affected 
by the gravitational forces from the star and from the disk. 
The planet's mass is set to 10M^ and no accretion onto the 
core is allowed. The core is placed at the corner of eight ad- 
jacent grid zones which are not treated specially in any way, 
so that the core itself is effectively unresolved by the grid. 
We assume that no circumplanetary envelope exists around 
the core in the initial state, so that initial state is defined by 
the background flow, unperturbed by the planet. In order to 
allow the later flow to adapt more gracefully to this condi- 
tion, we begin with a core mass of lJVf® and allow it to grow 
linearly with time to its final mass over the course of the 
first 10 years of evolution. We have not observed any side 
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effects of this adjustment over the rest of the simulation's 
lifetime. 

The simulation cube is centered on the core and is co- 
moving with it. The dimensions of the cube are set to ±47?h 
in each direction, corresponding to a cube of size ~ 0.897 AU 
on each side. In order to explore sensitivity to the core mass, 
we have performed one simulation with the core mass set to 
2OM0 , and we have performed several others with the cube 
volume as large as ±6i?H, depending on the goals of the 
simulation. Table [1] below defines the initial parameters for 
each simulation. 

The value of the x and z velocities (v x and fj z -each set 
initially to zero), corresponding to the radial and vertical 
velocities of the global initial condition, and the internal 
energy of the gas are mapped directly from the global initial 
condition onto the local condition. Since the grid spacing of 
the ID global initial state is preset to be identical to that 
on the corresponding 3D local grid, no interpolations are 
necessary. The values for the y velocity and volume density 
p require additional specification to complete the mapping 
from global to local. 

The local y velocity is set from the global azimuth ve- 
locities of the gas in the disk and the planet core. Because we 
assume a simulation cube that is fixed relative to the core, 
in order to obtain the y velocity, we must first subtract off 
the angular velocity of the frame 

v y = V4, — (dpi + a;)f2^,pi (21) 

where V$ is the orbit velocity of disk matter at each orbit 
radius and f2^ p i is the orbit velocity of the planet (and the 
rotation rate of the reference frame), and x is the offset from 
the planet's radial position, a p i. 

The volume density remains to be determined. In a com- 
plete model of a circumstellar disk, we expect to find gra- 
dients of most of the hydrodynamic quantities in both the 
radial direction and in the z direction, due to the near Kep- 
lerian character of the disk and to the radiative heating and 
cooling processes acting on it. In order to simplify our study 
however, we have chosen to neglect the z component of both 
stellar gravity and global disk self gravity, so that effectively 
no vertical structure in the disk exists. The symmetry is bro- 
ken only by the planet's gravitational force so that the flow 
around the core remains fully three dimensional. 

Since the global initial condition is defined by a surface 
density rather than a volume density, we require a conver- 
sion to completely specify the local condition. For this con- 
version we introduce the small inconsistency in our models 
that we use an estimate of the disk scale height (which of 
course depends on the existence of the z component of stel- 
lar and disk gravity), and the relation p — T,/H to specify 
the correspondence. The scale height is in turn defined by 
H — c 3 /Q and Q is the local rotation velocity in the disk. 
Its value of ~ 0.4 AU at the assumed orbit radius of 5.2 AU 
is approximately identical to the distance from the core to 
the edge of the simulation cube. 

2.3 Boundary Conditions 

In the shearing box formalism, the principle underlying the 
choice of boundary conditions is that no location in a disk 
is different than any other, excepting only that a shear term 
due to the background Keplerian motion must be accounted 



for explicitly. The box itself therefore represents one of many, 
essentially identical, replicas of a small subvolume of that 
disk. In this context, the common choice of implementing 
periodic boundary conditions for the shearing box (see e.g. 
lHawlev. Gammie fe Balbuslll995l ) means that a simulation 
of one small subvolume will be sufficient to represent the 
evolution of the entire disk. Clearly this principle does not 
hold for the present case, where our simlation volume con- 
tains a core and is therefore unique. 

Instead, we implement a slighly more restrictive princi- 
ple. Namely, we assume that the underlying disk background 
state is unchanging and is not affected by changes inside our 
simulation volume. We therefore set values for each state 
quantity in the boundary zones of the x and y coordinate 
directions of the local model (corresponding to the radial 
and azimuthal directions in the disk) to the initial values 
determined from the global model in the same manner that 
they would be had they instead been included in the simu- 
lation volume itself. 

For the y boundaries, the conditions are fixed to their 
initial values for the duration of the simulation. Because the 
simulation volume is assumed to be a cutout taken from a 
circumstellar disk, the y boundary quantities are not iden- 
tical at different x positions in the simulation, since they 
include the radial gradients of density, temperature and the 
rotation curve. This is an important consideration because 
at different x locations on the same boundary surface of the 
simulation box, corresponding to different orbit radii in the 
underlying circumstellar disk, y velocities may be directed 
either into or out of the simulation volume, with differing 
magnitudes. This behavior is a consequence of the initial 
near Keplerian flow of the underlying disk model, for which 
the y velocity varies as a function of the x position relative 
to the planet. 

Boundary conditions in the x coordinate are set assum- 
ing the same fixed initial condition. Flow into the grid is 
permitted but the value of each hydrodynamic quantity (in- 
cluding v x , the velocity perpendicular to the boundary) is 
reset to its initial value at each time step, so that the as- 
sumption of a fixed background state is maintained. Flow 
out of the grid is also allowed, should conditions require, by 
allowing v x to take on either zero or the value just inside 
the boundary if that velocity is directed outwards. 

This condition compensates for the development of 
large scale 'wakes', in the flow, corresponding to spiral struc- 
tures in the underlying circumstellar disk, which include 
some radial motion directed away from the planet in each 
direction. If a completely fixed condition is implemented, 
such structures can be reflected strongly and unphysically 
back into the simulation volume. Some reflection still occurs 
with the limited outflow condition, although at the level it 
is present, it proves essentially harmless because the amount 
of reflected material is both small and is quickly carried out 
of the grid through the y boundary. 

Initial simulations of mass propagating through simu- 
lation cube not containing a planet core demonstrated that 
inflow/outflow boundaries in the z coordinate resulted in 
large scale eddies forming in the flow. In order to avoid 
this problem, we use a reflecting condition on the z bound- 
aries instead. Employing this condition together with the 
ones for the x and y directions resulted in a quiescent flow 
throughout the simulation volume for as long (several hun- 
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dred years) as we cared to follow the evolution, and we did 
not pursue the origin of the eddies further. 



3 RESULTS OF SIMULATIONS 

Using the initial conditions and physical models described 
above we have performed a number of numerical experi- 
ments of the flow around a 10M^ core. One series of simula- 
tions vary the character of the background flow, by varying 
the components of gravity due to the disk that are included 
in describing the flow and by varying the background power 
laws of density and temperature. A second series of simu- 
lations vary the thermodynamic treatment of the gas, by 
changing the equation of state assumed for the gas and by 
changing the absolute scale of the background temperature. 
In table [T] we summarize several important parameters of 
the simulations. 

Columns 1-3 in the table specify a distinct label for each 
simulation and its resolution both in terms of the physical 
extent of the simulation cube, the number of grid zones and 
depth of the nesting. Columns 4-7 define the power law ex- 
ponents, p and q, for the surface density and temperature 
used to define the initial conditions of the global disk sys- 
tem, the assumed background temperature and a reference 
to the assumed equation of state for each model. Columns 
8-10 specify whether or not gravitational forces arising from 
the disk components described in sections 12.2.11 and 12.2.21 
have been included in the interactions with each other and 
the planet core in each simulation. The last column is the 
duration of each simulation in years. 

In the discussion below, we will first describe the evolu- 
tion of a 'prototype' model tm20, that is typical of the results 
obtained in many of our simulations. Then we will describe 
the results of our experiments that varied the rotation curve 
of the underlying disk, followed by experiments that varied 
the equation of state and background temperature of the 
disk. 

3.1 Morphological evolution of our prototype 
model 

In the discussion that follows and throughout this paper, 
we designate the model labeled tm20 as our 'prototype' 
model. Its features are typical of most of our simulations 
and the physical model underlying it is the most inclusive. 
This model was run with an ideal gas equation of state, with 
7 = 1.4 and included gravity from both the star and the core, 
as well as from both the global and local disk material. 

In figures [2}{5l we show 2D slices of the simulation vol- 
ume, taken through the disk midplane at a time 74 yr after 
the beginning of the simulation. In successive figures, we 
show the gas density, its entropy as defined by the expres- 
sion 

logS = logW, (22) 

its temperature and its velocity vectors projected onto the 
disk midplane. 

The density structures are highly inhomogeneous and 
become progressively more so closer to the core. Densities 
both above and below that of the background flow exist due 



to the development of hot 'bubbles' in the flow near the 
core, which then expand outwards. One such bubble is par- 
ticularly visible in the plots of the temperature distribution, 
emerging to the lower right. Such structures are common 
over the entire the duration of the simulation and emerge in 
all directions, depending on details of the flow at each time. 
Those details vary with time because of the very strong feed- 
back cycle that is generated as outgoing bubbles perturb the 
incoming gas flow that is responsible for generating later 
activity. Activity persists for as long as we have simulated 
the evolution without significant decay or growth. Also, al- 
though we have not included figures illustrating so explicitly, 
we note that plots of slices through our simulations in direc- 
tions perpendicular to the midplane show similar features: 
density structures which are highly inhomogenous, which 
become progressively more so closer to the core. 

Several ring-like structures surrounding the core are vis- 
ible and are remnants of intermittent large scale mass out- 
flows ('eruptions') from the environment of the core. These 
rings are sheared into an ovoid shape in the disk midplane, 
but are near spherical in the radial plane. The eruptions 
themselves appear at least in part to be the result of hot bub- 
bles that develop, expand and escape into the background 
flow and are due to shocks and compressional heating very 
close to the core. 

The structures within 1-2J?a of the core show little 
evidence of spiral arms seen in some pr evious work (e.g. 
iLubow et ai\\ 19991 ; iDAngelo et aZ.ll2003br i. but evidence for 
the development of other irregular dynamical structures is 
ubiquitous. Material orbiting the core or falling towards it 
encounters other material on different trajectories, generates 
shocks and heats the gas, causing bubbles and other com- 
plex morphological structures to develop. In actuality, such 
'orbits' are more often little more than single passes past 
the core as mass from the surrounding disk first falls into 
the core's gravitational influence, then returns to the disk 
and is swept away. 

We attribute the origin of the bubbles to the combina- 
tion of four conditions of our simulations. First, we include 
only a very small gravitational softening coefficient for the 
core, which means that the depth of its gravitational po- 
tential well is very deep and that the energy available for 
conversion into heat is very large. Second, we include no ex- 
plicit viscosity in our formulation of the fluid equations, so 
that features of the flow are not blurred or otherwise damped 
out by dissipative effects beyond our ability to model with 
precision. Also, our numerical method exhibits both very 
low numerical dissipation and very high fidelity for captur- 
ing the physical dissipation due to shocks that develop in 
the flow. Third, we do not include any mechanisms for en- 
ergy to be lost from the gas once it is converted to thermal 
energy from gravitational potential energy. The most im- 
portant such mechanism is of course, the losses that occur 
as the gas cools through radiation of photons. Finally, we 
do not include any mechanism for accreting material onto 
the core, because the dimensions of our finest numerical grid 
are comparable to the actual radius of the core. The ther- 
mal energy carried by that material is therefore also not 
accreted, remaining instead with the gas. In combination, 
these factors ensure thermal energy is produced efficiently, 
as flow features remain distinct enough to resolve the strong 
shocks that develop during the passage of material through 
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Local Box Around Proto Jupiter 

Resolution: 6x 64x 64x 64 Viewport Size= 0.897x 0.897x 0.897 AU 
Time Step: 31408 Time =7.456E+01 yr dt(l) =7.002E+04 s 
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Figure 2. The volume density in a 2D slice taken through the disk midplane for the full simulation volume of ±4i?H (top), and a blowup 
of the region within ±1/2Rjj OI the core (bottom) The spatial extent of each frame, in AU, is defined in its caption and labeled "Grid 
Size". Each cutout intersects the planet core at the origin. The white circles define the radius of the accretion sphere Ra = GM p i/c^ 
(small circle) and the Hill radius (large circle). The color scale is logarithmic and extends from ~ 10~ 10 to ~ 10 — 8 g cm -3 . 



© 0000 RAS, MNRAS 000, 000-000 



12 Andrew F. Nelson and Maximilian Ruffert 

Table 1. Simulation Parameters 



ijaDci 


Resolution 


Simulation 
Box Size 


P 


Q 






VjlODai 1J1SK 

on Planet 


^jlODai iJISK 

on Local Disk 


Local Disk 
on Local Disk 


T 

end 


lolO 


Q 
O 


X 


Ol 


±4/? H 


1 

1 


u 


1.0 


zUU 


k' 
IV 


eq.E 


no 


no 


no 


1600 yr 


ftio 





X 


D 1 


±4R H 


n 
u 





0.5 




IV 


eq.E 


no 


no 


no 


100 yr 


fllO 


e 
u 


X 


Ol 


±4R H 


n 
U 


r. 


1.0 


OOO 
ZUU 


IV 


eq.E 


no 


no 


no 


100 yr 


brlO 


rr 



X 


p. 1 3 
04 


±4R H 


i 


n 
U 


1.0 


zUU 


k" 
IV 


eq.E 


no 


no 


no 


100 yr 


bplO 


c 




X 


01 


±4R H 


i 


n 

U 


1.0 


900 


k' 

IV 


eq.E 


yes 


no 


no 


100 yr 


bslO 


5 


X 


64 3 


±4i? H 


i 





1.0 


200 


K 


eq.E 


yes 


yes 


no 


100 yr 


sglO 


5 


X 


64 3 


±4H H 


i 





1.0 


200 


K 


eq.E 


yes 


yes 


yes 


100 yr 


solO 


5 


X 


64 3 


±4i? H 


i 





1.0 


200 


K 


eq.0 


yes 


yes 


no 


100 yr 


sg20 


5 


X 


64 3 


±4i? H 


i 





1.0 


200 


K 


eq.E 


yes 


yes 


yes 


100 yr 


b05h 


6 


X 


64 3 


±6H H 


i 





1.0 


50 


K 


eq.E 


yes 


yes 


no 


100 yr 


tm05 


6 


X 


64 3 


±6i? H 


i 





1.0 


.50 


K 


eq.E 


yes 


yes 


yes 


100 yr 


tmlO 


6 


X 


64 3 


±5R H 


i 





1.0 


100 


K 


eq.E 


yes 


yes 


yes 


100 yr 


tm20 


6 


X 


64 3 


±4R H 


i 





1.0 


200 


K 


eq.E 


yes 


yes 


yes 


100 yr 


tm40 


6 


X 


64 3 


±4i? H 


i 





1.0 


400 


K 


eq.E 


yes 


yes 


yes 


100 yr 


isol 


6 


X 


64 3 


±4R H 


i 





1.0 


200 


K 


cq.HT] 


yes 


yes 


yes 


100 yr 


iso2 


6 


X 


64 3 


±4R h 


i 





1.0 


100 


K 


cq.[Tl] 


yes 


yes 


yes 


15 yr 


iso3 


6 


X 


64 3 


±4i? H 


i 





1.0 


50 


K 


cq.[[l] 


yes 


yes 


yes 


15 yr 


adil 


6 


X 


64 3 


±4fi H 


i 





1.0 


200 


K 


cq.\V2\ 


yes 


yes 


yes 


150 yr 



the core's environment, and dissipated (or simply lost, as 
in the case of accretion onto the core) inefficiently since no 
mechanisms are included to do so. Instead, the thermal en- 
ergy is converted again into kinetic energy as the high pres- 
sure gas expands outwards in bubbles and rejoins the back- 
ground flow. We will explore the consequences of energy loss 
mechanisms in section |4~41 below. 

Variations in entropy are significant in our simulations 
because the physical model includes only one mechanism by 
which its value can change after a packet of gas enters the 
simulation volume. Specifically, the hydrodynamic scheme 
itself generates entropy at shocks as a consequence of non- 
linearities in the solution of the fluid mechanics equations 
there. Therefore, any increase in entropy is, by assumption, 
a very strong a posteriori indication that shocks developed in 
the flow. Other entropy generation mechanisms that might 
be present in the flow will exhibit different characteristics 
than we observe in our simulations. For example, entropy 
generation in turbulent flows would be characterized by a 
conservative cascade of energy flow from large to small scale 
flow features, followed by entropy and thermal energy gen- 
eration at small scales. Such a cascade is not observed, and 
we believe turbulence is not an important characteristic of 
the flow activity present in our simulations. 

The fluid entropy in the background flow shows essen- 
tially no visible variation at different locations in the flow. 
Entropy of material near the core increases as it undergoes 
irreversible heating in shocks, but such heating is highly in- 
homogeneous and low entropy material is frequently present 
at distances only a small fraction of the accretion radius 
from the core itself. For example, comparing the bottom, 
high resolution, panels of figures [3] and [5] (discussed below), 
we see that low entropy 'background' material (i.e. the dark 
blue shaded material in figure^ in the upper right quadrant 
of the images is falling rapidly towards the core. Despite its 
proximity, it remains as essentially pristine background disk 
material until it falls to a distance of only < 0.1i?A (or equiv- 
alently, a few times Jupiter's current radius, Rj) from the 
core. There, its motion is strongly perturbed by the gravi- 



tational forces it experiences, which alter its trajectory onto 
paths that intersect those of other nearby packets of matter, 
generating shocks. As it re-enters the background flow, high 
entropy material begins to mix with the un-shocked back- 
ground material, but the mixed material remains distinct 
from the background until it leaves the simulation volume. 

Figure [4] shows the material temperatures for the same 
snapshots as seen in figurc[2] Temperatures as high as 6000 K 
are common in the innermost envelope, within distances of 
~ 0.1i?A(or as we note above, a few times Jupiter's current 
radius, Rj), but quickly drop to a few hundred degrees at 
larger separations. Most material at and outside the accre- 
tion radius, Ra, is only slightly warmer than the background 
200 K flow. Exceptions are materials that are clearly rem- 
nants of dynamical activity nearer to the core which expand 
outwards, such as the large mushroom shaped feature ex- 
tending to the lower right in the figure. Apart from the ra- 
dial temperature variation assumed for the background flow, 
visible temperature fluctuations of large magnitude do not 
extend away from the core to nearly the distances that are 
seen for the density variations. Nevertheless, material that 
has undergone a close encounter with the core does remain 
warmer and can be distinguished from the background flow 
even as it leaves the simulation volume, some time after the 
interaction occurred. 

It is of some interest to compare the distributions of 
temperature and entropy near the core. Although large de- 
viations are present, the temperature distribution appears 
far more spherically symmetric than does the entropy distri- 
bution. This is important because it indicates that material 
may undergo adiabatic compression without being shocked 
while falling inwards towards the core, and only suffer a 
shock deep in the envelope where the activity is most dy- 
namic. Implications of such compression and shock behavior 
will be discussed in section [S] 

Figure [5] shows velocity vectors of the flow in the mid- 
plane at the same time as is shown in figure [2] above. In con- 
flict to the orbital predictions from pure celestial mechanics 
that the flow of material approaching the Hill volume will 
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Figure 3. As in Figure [2] but showing entropy. The background flow is characterized by a small radial entropy gradient. Entropy 
increases by a factor ~ 6 over this value in the highly shocked neighborhood of the core, then decreases, as it mixes in with the low 
entropy background material in post-encounter evolution. 
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Figure 4. As in Figure [2] but showing temperature. The color scale is logarithmic and extends from 180 to 6500 K. Temperatures as 
high as 3-5000 K are common very close to the core, with temperatures decreasing rapidly to the background ~ 200 K value at increasing 
distance. 
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Figure 5. Velocity vectors shown projected onto the disk midplane on the coarsest grid in the top panel, and on the fourth nested grid 
in the bottom panel. 
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turn around on a 'horseshoe' orbit, matter approaching the 
outer portion of the Hill volume can be nearly unaffected, of- 
ten passing completely through its outer extent with only a 
small deflection. At the time of this snapshot, the trajectory 
of material radially outward of the planet passes directly 
through the Hill volume nearly unaffected, while material 
radially inwards experiences larger trajectory perturbations 
even at distances well beyond the Hill radius. Although we 
do not include specific examples as figures here, we note that 
as the flow evolves over time, this pattern frequently reverses 
itself, exhibiting exactly the opposite behavior-material in- 
ward of the planet is perturbed only slightly, while material 
outwards is significantly affected to a distances of 2-37?h- 
We will quantify these variations in more detail in section 
ET231 below. 

In contrast with this intermittently active flow further 
away, material is always very strongly perturbed on the scale 
of the accretion radius, where large amplitude space and 
time varying activity develops. This too conflicts with the 
orbital mechanics picture, in which matter inside the Hill 
volume simply orbits the core. Pristine disk material can en- 
ter the Hill volume from the background flow and shocked, 
high entropy material can escape and rejoin the background 
flow. Changes in the flow pattern occur on time scales rang- 
ing from hours to years, with a typical encounter time inside 
the accretion radius of less than a month. Additional discus- 
sion of these variations is found in section [3.2.31 below. 

3.2 Quantitative metrics of the activity 

The figures and discussion in the last section demonstrate 
the fact that the environment around the core is dynami- 
cally active in a qualitative sense, but do not quantify any 
measure of that dynamical activity apart from its overall 
morphology. Here we attempt to characterize the activity 
using a more quantitative analysis. First, we discuss the dis- 
tributions of mass and temperature in the envelope as a 
function of distance from the core, then, as integrated totals 
of mass and angular momentum contained in the envelope 
at a given time. We will use two metrics to describe the lat- 
ter two quantities. Specifically, we tabulate the total mass 
enclosed by a sphere of given radius and the total angular 
momentum of that material, calculated using the core as the 
origin of coordinates. 

3.2.1 The distribution of the gas near the core, and its 
temperature 

Because of the extremely non- homogeneous morphology and 
dynamically active flow, little purpose is served by attempts 
to reduce the flow variables to spherically averaged quan- 
tities. Nevertheless it remains interesting to analyze the 
overall distribution and limits of various quantities as they 
change as a function of radius. Such analyses are useful both 
to illustrate differences from the flow in our simulations from 
ID models and, secondly, to quantify the extreme values that 
may be expected at each radius and their significance in the 
context of the overall evolution of the system. 

Figure [5] shows the density and temperature of each grid 
zone in our prototype model, each as a function of distance 
from the core. In both cases, the distributions are nearly 
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Figure 6. The volume density (top) and temperature (bottom) 
of the gas in our prototype model, plotted for each zone in the grid 
as a function of distance from the core, at time t = 74 yr after 
the beginning of the simulation. The vertical dotted and short 
dashed lines define the location of and R\, respectively, and 
the long dashed line defines the present day radius of Jupiter, and 
is included to provide scale. Solid lines define the slopes of r~ 1,/2 
and r~ x / 4 power laws. 

flat at large distances (i.e. >>J?h), varying only slightly at 
a given distance as a consequence of the overall radial gra- 
dients of density and temperature defining the background 
flow in the global disk model. A quantity of material at large 
distances from the core exists that is clearly distinct from 
the background, as a population of points with somewhat 
higher temperatures, but is not visible in the density plot. 
This signature is due to material which has been processed 
through the core's envelope and returned to the background 
flow, but which has not yet fully cooled or mixed with it. 
The material retains its higher temperatures until it is car- 
ried out of the grid volume. An example of such flow is illus- 
trated in figures 121 II where a hot outflow is shown emerging 
to the lower right in the lower panels of the figures. 

At distances smaller than the Hill radius, both density 
and temperature begin to deviate more significantly from 
their background values. At a given distance from the core, 
the maxima found in either quantity approach factors of 
~ 5 — 8 larger than the minima at the same distance, with 
the largest differences occurring inside the accretion radius. 
Densities both increase above their background values and 
decrease to well below the background, reflecting both the 
compression of infalling material and the expansion of hot 
material after it is processed deep in the core's gravitational 



© 0000 RAS, MNRAS 000, 000-000 



Dynamics of Core Accretion 17 



potential well. Both quantities increase all the way to the in- 
ner radial limit imposed by our grid resolution, at a distance 
from the core of ~ 5.6 x 10 9 cm, which is somewhat smaller 
than the present day radius of Jupiter itself. For the inner- 
most zones, densities reach values of ~ 10 -8 g/cm 3 , nearly 
two orders of magnitude above that in the background flow, 
while temperatures rise to 5-7000 K. Examination of lower 
resolution variants of this simulation (i.e. simulations sglO 
and lolO, not shown here), show lower maximum densities 
and temperatures, indicating that the values reached in fig- 
ure [6] are not converged, and suggesting that still higher 
values might be reached near the core/envelope interface if 
resolution were to be increased. 

Both distributions vary widely at any given radial dis- 
tance from the core, and neither can reasonably be fit to a 
ID power law. In particular, the density distribution can- 
not be fit to a r~ 2 power law, as would be expected from a 
spherically symmetric infall. Instead, and to the extent that 
it follows a single profile at all, it falls nearer to a profile 
between r _1//4 and r -1 / 2 , as indicated by the lines accom- 
panying the figure. The much shallower radial dependence 
presumably reflects the fact that no mass accretion onto the 
core is permitted, so that material entering the envelope is 
opposed by material already there, and by material return- 
ing to the disk. 

The lower limits of the temperature distribution at each 
radius provide an interesting flow metric because they de- 
fine the population of material at each radius which has 
undergone the least heating, either due to compression or 
to shocks. Of great interest is that material at tempera- 
tures below ~ 300 K is not uncommon even at distance as 
small as 10 11 cm, about 1/4 of the accretion radius, while at 
other locations at the same distance, temperatures as high 
as 1-2000 K are present. The presence of low temperature 
material close to the core demonstrates that some fraction 
of the background disk material can penetrate deeply into 
the envelope while undergoing essentially no processing by 
it, either due to shocks or to adiabatic compression. This 
material is interesting because as it continues its trajectory 
through the envelope, it may subsequently undergo some 
compressional or shock event which heats it rapidly to tem- 
peratures more typical of other material at these distances. 
A quantification of what fraction of material undergoes such 
rapid heating, as opposed to that undergoing some slower 
heating process or to that simply being advected out of the 
envelope again, is beyond the scope of this study. Its exis- 
tence is important to note however, because if an appreciable 
quantity of material can penetrate deeply into the envelope 
while remaining cold, then be rapidly heated and ejected 
again into the background flow, some signature of the pro- 
cessing may remain in material available for study today. We 
discuss a preliminary study of such a possibility in section 

El 



3.2.2 The mass contained in the envelope as a function of 
time 

Figure [7] shows the mass enclosed by spheres of radii, Ra 
and r = 0.10 AU~ 0.9i?H, each as functions of time. Except 
for the first ten year period of the simulation (omitted in the 
plots), in which we artificially grow the core to its final mass, 
the enclosed masses exhibit no secular mass accretion over 
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Figure 7. The mass enclosed by spheres of radius r = 0.10 AU~ 
0.9-Rjj (top) and r =Ra(= 0.026 AU) (bottom), each as functions 
of time during the simulation. The first 10 years of the evolution 
define the period when we first 'grow' the core from IMq to its 
final mass of lOAfgj , thus the envelope masses at this time are not 
meaningful and are not included in the figure. The dotted lines 
correspond to the enclosed mass at t = 0, as characterized by the 
unperturbed disk background condition. 



the duration of the simulation, consistent with our neglect of 
cooling. The amount of mass enclosed within the accretion 
radius is quite small, averaging ~ .007M® , while on the scale 
of the Hill radius, a quantity of ~ 0.24M® collects near 
the planet. In comparison to the 'background flow' initial 
condition, the inner envelope mass has increased by a factor 
~ 1.8, while the envelope mass as a whole has increased by 
a factor ~ 1.2. 

Variations in the enclosed masses with time are visible 
in both curves, with magnitude ~ 10 — 20% for the smaller 
sphere and ~ 5 — 10% for the larger. In neither case does the 
variation grow or shrink as a function of time, again consis- 
tent with our neglect of cooling. The temporal character of 
the variability, such that any sort of periodicity can be de- 
fined at all, clearly favors higher frequency oscillations close 
to the core and slower oscillations farther from it. Variations 
on time scales well under a year are common for the accre- 
tion sphere volume, and are often superimposed on longer 
time scale features, as the overall morphology of the flow 
reacts to changes in the larger scale flow. The existence of 
the variations over time is a consequence of the dynamical 
activity discussed above, as morphological features flow into 
or out of the respective spheres defining our metrics. In turn, 
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their time scales are consequences of the shorter dynamical 
times deeper in the core's gravitational potential well, as 
compared to those farther from the core. 

3.2.3 The angular momentum of the gas around the core 

In addition to the mass and temperature distributions dis- 
cussed in the last sections, the distribution of angular mo- 
mentum of the material in the environment of the core can 
also provide insights into the character of activity. Here, we 
describe the magnitude and the variability in angular mo- 
mentum distributions of material around the core. In the 
figures below and the discussions accompanying them, we 
will show both instantaneous values as functions of time, 
and time averages and RMS deviations over the duration of 
each simulation. The time dependent plots will provide in- 
sights into any periodicities that are present in the results. 
Time averaged quantities will be useful for comparisons be- 
tween simulations with different physical parameters, and 
in comparison to the present day values for Jupiter. To en- 
able such comparisons, we will present values normalized by 
envelope mass, as specific angular momenta. 

In the discussion below, we refer to the envelope's angu- 
lar momentum as its 'spin'. We caution the reader however, 
that this term is not strictly correct, since while some ma- 
terial may be in a bound orbit around the core, other parts 
of the 'envelope' for which spin is tabulated, may instead be 
more accurately classified as part of the background flow. 
Also, as an aid to reader's intuitive grasp of the meaning of 
each of the three spin components, we note that the z com- 
ponent corresponds to the 'top'-like motion of the planet, 
and corresponds to the largest component of the present day 
spin of the planets in our solar system (excluding Uranus). 
The x component defines a spin axis outwardly directed 
along the line connecting the star and the core and the y 
component defines a spin axis pointing forwards along the 
direction of motion of the core in its orbit. 

In figure [8] we show the enclosed spin angular momen- 
tum of gas relative to the core as a function of time, in 
each of the three coordinate directions, and at the same 
two radii as the enclosed mass above. For the larger sphere, 
the spin magnitude oscillates, seemingly aperiodically, both 
above and below zero in the the x and y directions on time 
scales of 10-20 yr. The peaks and valleys of the oscillations 
extend to values as large as a few times that of Jupiter's spin, 
normalized by envelope mass, and typically falls in a range 
within ±1 Sj of its mean. In the z coordinate, a similar fea- 
ture, with oscillations of similar or slightly larger magnitude, 
is also present but does not cause the sign to change because 
the mean lies significantly below zero. The fact that the sign 
of the spin is uniformly negative is important because it in- 
dicates that the 'envelope' contains material which is merely 
passing through the Hill volume as part of the background 
shear flow. It is not a part of either a static or forwardly- 
rotating envelope, as expected in a purely dynamical flow in 
which pressure plays no role. 

On the smaller scale of the accretion radius, where we 
would expect the matter flow to be influenced more by the 
core, the enclosed spin per unit mass again lies near par- 
ity with present day Jupiter. Its magnitude varies on time 
scales much smaller than a year, both above and below zero 
(i.e. both one direction and the other), in all three coordi- 



nate directions. The time variations for this smaller enclos- 
ing sphere do not appear to be well correlated with any of 
the longer variations visible in the larger sphere. In large 
part, this is a consequence of the much larger total mass 
enclosed by the latter, so that the variations on small scales 
are simply washed out. 

Combining the results of figure [7] which shows that the 
envelope contains only a tiny fraction of the present day 
mass of Jupiter, with those of figure [5] which shows spin 
normalized to the envelope mass, we find that actual spin 
magnitudes are only a small fraction of present day Jupiter's 
total spin. At such an early stage in the evolution, such a 
value is not particularly surprising. Close examination of 
the spin characteristics remains interesting nonetheless, for 
several reasons. First, perhaps the most important feature 
of figure [8] is not a exact values that any spin component 
takes, but rather that the amplitudes of each component 
(apart from the background shear of the disk included in the 
z spin calculation), are quite similar. This indicates that the 
spin direction of the planet itself is not yet fixed. The era 
of giant impacts between nascent cores may still continue 
during this phase of planet formation, with little evidence 
remaining in the spin characteristics of the planet in its final 
configuration. 

Secondly, the similarity between the spin magnitudes 
seen here and those of Jupiter indicates that no angular 
momentum loss mechanism must be invoked in the circum- 
planetary environment in order for accretion to proceed. 
The most common such mechanism would, of course, be a 
circumplanetary accretion disk through which circumstel- 
lar disk material would pass as it accretes onto the core. 
No such disk appears neccessary at this early stage of the 
planet's evolution. 

Third, the envelope's spin direction does not appear 
to be established: it changes sign in all three coordinate 
directions on spatial scales comparable to _Ra and it changes 
sign in x and y even on much larger scales. Also, although the 
z coordinate spin in the r = 0.1 AU enclosing sphere varies 
widely, it varies around a negative mean value. At this scale, 
a purely dynamical flow (i.e. that of a set of 'test particles', 
affected by gravitational and fictitious forces only) , would be 
expected to produce forward envelope rotation (positive s z ) 
as mass continues to accrete. Instead, some portions of the 
background material flow through the Hill volume, because 
the purely dynamical picture is incomplete. The flow is not 
dominated solely by the gravitational forces of the core and 
star combined with the fictitious forces present in a rotating 
medium, but also by the pressure forces of the gas itself. We 
will discuss the relative importance of the hydrodynamic 
effects in more detail in section 14.4.21 below. 



4 SENSITIVITY OF THE ACTIVITY TO 
NUMERICAL AND PHYSICAL 
PARAMETERS 

The dynamical activity shown in the results above persists as 
long as we are able to run our simulations-typically about 
100 years of simulation time for each run. Unfortunately, 
that time is quite short compared to the time scale for the 
envelope to grow in mass or for the disk as a whole to evolve. 
We are therefore unable to answer many important ques- 
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Figure 8. The spin angular momentum per unit mass calculated around the core for gas enclosed by a sphere of radius r = 0.1 AU(~ 
0.9-Rh) or by a sphere one accretion radius (r =Ra(~ 0.026 AU) in size (top and bottom respectively). For the larger sphere, the net 
2 spin has a negative bias due to the shear of the background flow. Both the x and y spins are centered near zero. All three exhibit 
large amplitude swings in magnitude, including changes in sign: the sense of rotation around the core reverses itself. In contrast, the spin 
components enclosed by the accretion radius exhibit almost no systematic trends except a continuing large amplitude and short term 
variation. 



tions regarding the activity. In particular, in order to make 
broadly applicable conclusions about the physical systems 
we are modeling, we must show that they are common to a 
wide variety of initial conditions that may be encountered 
in the circumstellar nebula at various times in its evolution, 
that they are common to a variety of physical models, and 
that they are not unduly contaminated by consequences of 
our initial conditions. The results must not simply be 'start 
up transients' or other artifacts of initial conditions that do 
not represent the true system with sufficient fidelity. 

To fully address these concerns requires models far more 
sophisticated than are presently possible and we defer dis- 
cussion of any final determinations of their answers to the 
future. Instead, in this section we consider a number of sim- 
ple parameters that can be varied within the context of our 
models, and to which we may reasonably suspect that our 
models may be sensitive. If our models vary their behavior in 
these, admittedly crude, sorts of parameter studies, we may 
justifiably conclude that they will be even more sensitive 



to the much more varied and variable physical properties 
present over the lifetime of a circumstellar disk. In addi- 
tion, there may be unaccounted for numerical issues beset- 
ting our simulations which cast doubt on the results. If the 
simulations are comparatively insensitive to the parameters 
varied here, we will make a first step towards the physically 
meaningful conclusion that the dynamical activity discussed 
above can be present at interesting levels at many different 
times during circumstellar disk evolution. 

As a useful means to quantify similar and dissimilar 
behavior, we adopt the spin of the envelopes, introduced in 
the last section, as a metric suitable for our purposes here. 
We summarize the time averaged quantities for each of the 
simulations discussed above and in the following sections, in 
table(2] As in table[T] the first column defines the simulation 
name. The next two pairs of columns (columns 2/3 and 4/5) 
each specify the time averaged spins and their variances of 
the envelope material enclosed by spheres of r = 0.1 AU 
and r =Ra, respectively, for each simulation. Both quan- 
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Table 2. Time Averaged z components of Spins and Spin Vari- 
ances 



Label 


Sz/sj 


Ca 


Sz/sj 


fa 


CR offset 




r = 0.1AU 


r = 0.1 AU 


r =R A 


r =R A 


5a (Ru) 


lolO 


-3.54 


0.99 


0.28 


0.73 


-0.219 


ft 10 


-4.63 


1.31 


-0.05 


1.20 


-0.039 


fllO 


-4.40 


0.91 


0.005 


0.97 


-0.172 


brlO 


-4.75 


0.93 


-0.43 


1.19 


-0.203 


bplO 


-4.83 


1.31 


-0.53 


1.15 


-0.344 


bslO 


-4.48 


1.27 


-0.13 


1.19 


-0.203 


so-1 D 
SglU 


-4 68 


0.98 


-0.18 


1.13 


-0.203 


solO 


-3.99 


0.07 


-0.17 


0.03 


-0.203 


sg20 


-3.36 


1.71 


-0.03 


1.93 


-0.156 


b05h 


-2.75 


1.70 


0.60 


1.67 


-0.053 


tm05 


-2.26 


1.50 


0.78 


1.54 


-0.053 


tmlO 


-4.26 


1.43 


0.11 


1.41 


-0.103 


tm20 


-4.37 


1.24 


0.038 


1.05 


-0.203 


tm40 


-5.53 


0.76 


-0.52 


0.82 


-0.39 


isol 


-2.59 


0.31 


1.44 


0.42 


-0.203 


iso2 










-0.102 


iso3 










-0.051 


adil 


-3.88 


0.13 


-0.22 


0.15 


-0.203 



tities normalized to the angular momentum per unit mass 
of present day Jupiter. For the series of simulations with 
varying background temperature (tm05-tm40), the sphere 
of radius _Ra, is replaced by one of radius r = 0.025 AU, for 
consistency, as discussed in section 14.4.21 For the iso2 and 
isoS simulations, no time averages are reported, due to their 
short duration. The last column specifies the radial offset of 
the corotation radius from the planet that results from each 
physical model. We explore the importance of variations in 
this quantity below. 



4.1 Sensitivity to the changes in the rotation 
curve 

The dynamical activity seen above is clearly driven by inter- 
actions between the background flow and the core, mediated 
by hydrodynamic interactions deep in the core's gravita- 
tional potential well. We saw that the background flow could 
be distinguished from the local envelope to much smaller dis- 
tances than are predicted by purely dynamical models. We 
also saw that background material could penetrate deeply 
into the core's potential well before undergoing substan- 
tial thermodynamic evolution such as compression or shock 
heating. It is therefore reasonable to ask whether features 
of the background flow may influence the strength of other 
properties of the activity. Does changing properties of the 
background flow change the behavior of activity generated 
in our simulations? 

Two important properties of the background flow are 
first, the overall shear itself and, second, the radial location 
of 'corotation' relative to the orbit radius of the core itself. 
The first is largely a consequence of the near Keplerian or- 
bital flow in the disk and, while other rotation curves may 
be of interest in some systems (e.g. galaxies), they are of 
no practical importance in the present context. The relative 
positions of the core and of the disk material orbiting at 
the same frequency is of much greater interest. It may shift 
signficantly due to the overall mass of the disk (and whether 



that mass is included self consistently in a calculation of the 
rotation curve), its temperature profile and its density pro- 
file. To what extent will such shifts affect the flow around 
the planet itself? 

Although only a model that includes all components 
of gravity on all components of the system can be correct, 
including, excluding or otherwise modifying one term or an- 
other will allow us to explore the sensitivity of the results 
to that term. Here, we will include the effect of gravity and 
pressure in several different ways in order to explore the 
sensitivity of the results to the importance of the rotation 
curve of the disk, and the offset it may have relative to 
the orbit velocity of the planet. All experiments will include 
the gravitational force of the star and the planet as a fixed 
background. In separate simulations we will add additional 
sources of gravitational force, component by component. Pa- 
rameters for each of these simulations, and the position of 
the co-rotation resonance relative to the core for each of 
these simulations are tabulated in Tables [T] and [2] 

In figures [9] and \W\ we show the z component of spin 
per unit mass enclosed inside spheres of the same two radii 
around the planet as are shown in figure [8] for a variety of 
different physical assumptions about the underlying rotation 
curve. In every case, variations qualitatively similar in both 
magnitude and time scale to those seen for our prototype 
model are present. Quantified in terms of their variances 
over time, the spins in each of the models vary by as much 
or more than the present day normalized spin of Jupiter over 
time scales much shorter than the 100 yr lifetime of our sim- 
ulations. As expected from their different initial conditions, 
differences appear in the time averaged spin values, with a 
spread of about half of the present day normalized spin of 
Jupiter between the maximum and minimum seen in a given 
model. Of note is that the spread in the average values is 
smaller than the variations seen over time in a single model. 
No model stands out as particularly different from the oth- 
ers, nor are the results of any model particularly different 
from that of the prototype model discussed in section 13.11 
In fact, the difference between the prototype's averages and 
those of the sglO simulation, which uses an identical physical 
model but one less grid, are comparable to the differences 
between it and any of the other models with varying physical 
parameters. 

We conclude that different physical conditions in the 
disk will generate difference in the details of the envelope 
activity, but will not affect its presence or overall magnitude. 

4.2 Sensitivity to the description of the core itself 

In addition to possible sensitivities to the background flow, 
we might also expect dynamical activity to be sensitive to 
the mass of the core and to the distribution of mass close to 
the core, each through the gravitational forces they exert. 
Both can readily be probed by varying the mass of the core 
and its gravitational softening length. Increasing the core 
mass retains the shape of the potential well, but changes its 
depth. Increasing the softening length changes the shape of 
the potential well, while decreasing its depth (see eq. I14[) . 
The softening length serves as a probe of the sensitivity to 
the mass distribution because, in addition to its use to re- 
duce numerical instabilities in a simulation, its presence in 
the otherwise perfectly inverse square, Newtonian force law 
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Figure 9. The z spin angular momentum per unit mass calculated around the core for gas enclosed by a sphere of radius, r = 0.1 AU(~ 
0.9-Rh) in size, for the model designated in the upper left corners of each panel. The time averaged magnitude of the spin and its variance 
are tabulated in the upper right corner of each panel. 
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Figure 10. Same as for figure llOl but for a sphere of radius r =-Ra(~ 0.026 AU). 
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is effectively equivalent to the statement that the core's mass 
distribution is spatially extended on the scale of the soften- 
ing length itself. 

Here, we test the sensitivity of our results with two mod- 
els. In one model we define a core mass of 20Me , which is 
twice our prototype value, and in the second model, we de- 
fine the gravitational softening parameter, e feq. I14[l . from 
e = 1 x 8x to e = 16 x Sx, where 8x is the size of one grid 
cell on the finest (in this case, the fifth) mesh. This increase 
corresponds to an increased spatial scale from ~ 1.9Rj, to 
~ 30Rj, or ~_Rh/8. Figure [TT1 shows the z spin of envelope 
material for these two models. In both cases, the character- 
istics of the envelope activity differ markedly from those of 
the prototype. With increased softening, the variations are 
near completely absent. With increased core mass, the mag- 
nitude of the variations increase by a factor of available to 
disk material as it traverses through the deeper gravitational 
potential well of the more massive core. 

We conclude that characteristics of the core and en- 
velope do affect the magnitude of the dynamical activity, 
and that a deep, sharp potential well is a requirement for 
the activity to be present. An envelope will likely exhibit 
dynamical activity only when the core it surrounds is suffi- 
ciently massive, and when the envelope itself is sufficiently 
diffuse not to contribute substantially to the shape of the 
potential well. 

4.3 Durability of the dynamical activity 

Because the duration of our simulations is short and the ini- 
tial conditions are not especially tuned to begin close to any 
sort of equilibrium or steady state flow, it is possible that the 
dynamical activity is simply a consequence of inconsisten- 
cies in our initial condition. If so, then we might expect the 
behavior of the flow observed in a given simulation to change 
as a function of how long we run that simulation. In order 
to investigate this possibility, we have run a comparatively 
low resolution simulation, using the same physical model as 
our prototype, for more than 1600 years of simulation time. 

In figure [12] we show the z spin of the envelope mate- 
rial as a function of time, over the life of the simulation. 
As for the high resolution models above, the magnitude of 
the spin displays continuing activity, with no signs of de- 
crease or modulation. Inside the accretion radius, the spin 
direction changes on time scales of < 1 yr, as it does for 
our prototype. On the scale of the Hill volume, the over- 
all direction (backwards relative to the present day spin of 
Jupiter) again retains a strong signature of the background 
flow, but is again perturbed significantly by the dynamical 
activity. In both cases, the variances around the means fall 
some ~ 20% smaller than were seen in the higher resolution 
models, but remain comparable to the normalized spin of 
present day Jupiter. An important difference is that a posi- 
tive mean spin exists for the material inside _Ra and, for the 
r < 0.1 AU volume, a somewhat less negative overall spin 
because of the opposing contributions of the inner envelope 
and the background flowing through the outer envelope vol- 
ume. 

We attribute the differences in behavior to differences 
in resolution because the smaller cells and time steps pro- 
vided by the higher resolution realization permit more rapid 
and larger amplitude variations than in the lower resolution 



simulation. None of the differences contribute to any over- 
all growth or decay of the dynamical activity. We therefore 
conclude that the dynamical activity is not due to start up 
transients arising from inconsistencies in our initial condi- 
tion. 

4.4 Sensitivity to thermodynamics 

We now turn to the question of how our results depend 
on the thermodynamic treatment of the gas. There are two 
distinct physical properties that may be important for the 
flow. They are first the heating and cooling properties of 
the gas together with the equation of state and second, the 
temperature of the background flow. In the following two 
sections, we examine both of these possibilities in turn. 

4-4-1 Importance of heating and cooling, as approximated 
by assumed equations of state 

Although we do not include any external heating or cooling 
directly in our study, we can perform experiments compar- 
ing our ideal gas equation of state to locally isothermal and 
locally isentr o pic eq uatio ns of state. Detailed discussions in 
iNelson et all l|200d ) and IPickett et all l|2003h describe the 
physical interpretation of the effective heating and cooling 
that can be placed on simulations that use such formula- 
tions. For our purposes here, it is sufficient to note that both 
ordinarily imply very efficient cooling. Simulations which 
employ a locally isentropic equation of state effectively as- 
sume that all shock heating is instantly radiated away. Sim- 
ulations which employ a locally isothermal equation of state 
effectively assume that all heating due both to shocks and 
PdV work is instantly radiated away. 

Figure [131 shows the density structure obtained from lo- 
cally isentropic and isothermal evolution. In each case, the 
morphologies are far more uniform than are seen in the evo- 
lution with an ideal gas treatment (as seen in figure [2}. The 
density structure in the locally isentropic evolution appears 
nearly spherically symmetric, with only small perturbations 
being visible in the image, towards the upper left and lower 
right of the core. Over time, perturbations of similar size 
and magnitude appear at various positions in the envelope 
as the flow changes, but do not grow substantially larger. 

Locally isothermal evolution leads again to fewer dy- 
namical features in the flow than in our prototype simula- 
tion, but here some activity remains visible, predominantly 
in the form of a barred spiral pattern generated by material 
orbiting the core, along with some smaller features which 
are remnants of prior spiral structures that have merged or 
disintegrated as the evolution proceeds. 

The densities in the deepest parts of the envelope in 
both of the simulations rise to values much higher than those 
seen in our prototype. The highest central densities occur in 
the isothermal evolution, at more than four orders of magni- 
tude higher than in the ideal gas evolution, and two orders of 
magnitude higher than the isentropic evolution. The higher 
densities are the most visible consequences of the implied 
cooling assumptions in each model, which do not permit 
high temperatures to evolve deep in the envelope. In conse- 
quence, pressure forces opposing the gravitational attraction 
of the core are correspondingly reduced, which allows addi- 
tional mass to concentrate there. 
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Figure 11. The z spin angular momentum per unit mass of the envelope material enclosed in spheres of radius r = 0.1 AU(~ 0.9-Rh) 
(top) and of r =i?A (bottom), each as functions of time. The left panels refer to simulation so 10, with increased gravitational softening, 
while the right panels refer to simulation sg20, with increased core mass. 



We monitored the stability of the gas to numerically 
induced collapse thr ough violation of the Jeans criterion 
l|Truelove et aZJIl997f ) and, although not strictly applicable 
to this systenQ the criterion was satisfied over the evolu- 
tion of the simulation, using a safety factor of J ~ 4. Nev- 
ertheless, at time t w 85 yr, the gas near the core began to 
collapse in upon itself, and shortly thereafter we terminated 
the simulation. Although a strict application of the Jeans 
criterion is satisfied, we believe that this collapse behavior 
is not due to any physical cause, but rather a numerically in- 
duced consequence of the isothermal evolution and the high 
spatial resolution. 

Figure [14] shows the evolution of the z spin of the en- 
velope material over time, for both simulations. In both, 
some variability remains, particularly in the isothermally 
evolved model. In comparison to the variations seen in fig- 



The Jeans analysis considers small, linearized perturbations of 
quantities away from a smooth background, while the present 
configuration also includes a large 'external' potential gradient, 
due to the core. 



ure [8] however, the magnitudes are far smaller. Variations 
of only <0.2 Sj are found in the isentropic evolution com- 
pared to ~ 1 Sj or more seen in figure [8] The time scales 
of the variations also appear somewhat longer, even for the 
smaller volume accretion sphere, defined by r < Ra, though 
we have made no attempt to quantify such time dependence. 

Variations seen in the isothermally evolved simulation 
are larger than the isentropic simulation, near ~ 0.3 — 0.4SJ, 
but again remain smaller than those seen in the prototype. 
For the larger volume sphere, the spin remains consistently 
negative, as it does for all of our other simulations, but its 
magnitude is much smaller. This difference appears to be 
due to the fact that the spin of the inner envelope is consis- 
tently positive, acting in the opposite sense to the flow fur- 
ther from the core. Quantified by its behavior for r <Ra, the 
sign of the envelope spin is consistently positive: though its 
magnitude still varies, the spin of the envelope is always ori- 
ented in the same direction as the present day Jovian plan- 
ets. Positive spin orientation is an expected result for purely 
dynamical systems of course (i.e. those without hydrody- 
namic effects), because of the influence of the fictitious forces 
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Figure 12. The z spin angular momentum of the envelope material in simulation lolO as a function of time, enclosed in spheres of radius 
r = 0.1 AU(~ 0.9-Rh) (top) and of r =Ra (bottom). The left panels show the variation over the entire span of the simulation, while the 
right panels show the variation over a 100 yr segment. 



on the motion. The fact that we recover a positive spin ori- 
entation only in this simulation, where hydrodynamic feed- 
back is largely suppressed through the isothermal evolution 
assumption, provides strong support for our conclusion that 
the presence and strength of the dynamical activity is tied 
directly to hydrodynamic feedback effects on the flow. Fur- 
ther, it serves as an important check on the physical model 
underlying all of our simulations because under physical con- 
ditions where feedback is removed, the expected dynamical 
result is recovered. 

Due to the changes in the behavior of the flow in these 
models, we conclude that assumptions made regarding the 
energy balance will be of critical importance for models 
of dynamical activity in the core's environment. Of these 
assumptions, a complete model for the radiative transport 
through the gas will be of primary importance deep in the 
envelope. A necessary corollary is that a well determined 
dscription of the material composition and size distribution 
over the range of temperatures and densities will also be crit- 
ical, because of their influence on the opacities throughout 
the relevant envelope volume. 



4-4-2 Importance of the background temperature 



In section 13.2.31 we concluded that an important physical 
parameter in defining properties of the activity was the 
existence of hydrodynamic feedback in the system, arising 
from pressure forces. As the pressure and pressure gradients 
change in importance compared to other forces acting on 
the flow, will the character of the activity also change? Here 
we explore this question. 

The most physically relevant parameter to vary in such 
a study is of course the pressure itself, or its gradient, but 
neither are directly accessible for use as numerical param- 
eters in our simulations. Instead, we vary the overall back- 
ground temperature assumed for the disk in different sim- 
ulations. This quantity will serve as a proxy for the overall 
magnitude of pressure effects through the relationship de- 
fined by the equation of state. To address our question, we 
have run in four variants of our prototype model (simula- 
tions tm05, tmlO, tm20 and tm40), with background tem- 
peratures at the core's orbit radius of T = 50, T = 100, 
T = 200 and T = 400 K. For the T = 50 K model, we 
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Figure 13. As in FigureO but showing densities for the locally isentropic (top) and isothermal (bottom) simulations. The color scales 
are logarithmic and extend from ~ 7 X 10 — 11 g/cm 3 in each panel to ~ 10 -7 and ~ 2 X 10 -5 g/cm 3 for the isentropic and isothermal 
simulations, respectively. The black circles define the accretion radius, for each simulation. 
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Figure 14. The z component of spin angular momentum of the envelope material, for the same two spherical volumes shown in figure 
[8] and for the isentropic (left panels) and isothermal (right panels) equation of state assumptions. Evolution of the isothermally evolved 
simulation suffered catastrophic numerical errors after ~ 80 yr of evolution and was terminated (see text). Average values appearing in 
each panel are derived from the period spanned by [10,80] yr for the isothermally evolved simulation and [10,150] yr for the adiabatically 
evolved model. 



also include a variant omitting local self gravity of the disk 
(b05h). 

We include two variants of the model with T = 50 K 
background temperature in order to investigate the conse- 
quences of a numerical issue affecting simulation tm05. Due 
to its low background temperature and the inclusion of lo- 
cal self gravity, mass accumulates to an unphysical extent 
inside the simulation volume. We attribute this behavior to 
the fact that the Jeans wavelength derived for the material 
within the simulation volume appoximates the size of the 
simulation volume itself, rather than to any physical cause. 
Though an important consideration for other metrics, such 
as enclosed mass, in terms of the effects this numerical de- 
fect may have on the metric chosen in the present discussion 
(spin), we will find that the differences between the two sim- 
ulations are minimal. 

Figure [15] shows the evolution of spin of the envelope 
material in the z coordinate over time, for each of these sim- 
ulations. From left to right, the figure shows the progression 
from low to high background temperatures. For consistency 
with the size of the sphere defined for the accretion radius 



in previous sections, we plot the spin enclosed by a sphere 
of radius r — 0.025 AU, rather than r =Ra, since the lat- 
ter quantity varies with temperature. The former quantity is 
approximately the value of Ra at the T = 200 K background 
temperature of our prototype simulation. 

Several notable features appear in the figure. First, the 
time averaged spin values vary with background temper- 
ature. For the smaller enclosing sphere in particular, the 
highest background temperature models generate spins av- 
eraging well below zero (opposite the present day spin di- 
rections of the major planets). For the simulation with the 
lowest background temperature, the spin increases to an av- 
erage nearly 80% that of Jupiter's present day spin normal- 
ized by mass and in the same direction. At larger radii, a 
similar effect occurs-the envelope spin increases, becoming 
much less negative as temperature decreases. The values of 
the spin averages fall well outside the range derived from the 
study of different physical models in figures [9] and [TO] above, 
and we conclude that they signify a quantitative difference 
in the flow behavior in comparison to the other models in 
our study. 
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Figure 15. The z component of spin angular momentum per unit mass as a function of time, enclosed in spheres of radius r = 0.1 AU(~ 
0.9-Rjj) (top) and of r = 0.025 AU (bottom). From left to right, the panels show the simulations with background temperature T = 50 K, 
T = 100 K, T = 200 K and T = 400 K, respectively. In the left-most panels, simulation tm05 is shown with solid curves, and b05h is 
shown with dotted curves. 



The trend towards more negative spins at higher tem- 
peratures and more positive values at lower temperatures 
correlates directly with the changes in position of co-rotation 
relative to the core (see Table [2}. For simulation tm40, for 
example, the inwards offset is largest of any in our study. By 
its nature, such an offset exposes the core to material with 
correspondingly larger relative velocities, and larger net an- 
gular momenta relative to the core. In this case, the large 
negative angular momentum of the material overwhelms the 
tendency of the system's fictitious forces to produce pro- 
grade rotation around the core, even in the regions closest 
to the core itself. 

At the other extreme, the behavior seen in the lowest 
temperature simulations is consistent with the prograde ro- 
tation we would expect in the limit of a zero temperature, 
purely dynamical flow, expected from an ensemble of mass- 
less test particles. Even so, for the T = 50 K simulation, 
where Ra^Ru, the spin enclosed within the Hill volume re- 
mains negative, indicating that a significant component of 
material at these distances is not 'envelope' material at all, 
but rather merely background disk material which is passing 
through that volume of space. We conclude that for any tem- 



perature that we might reasonably expect to be present in 
the disks of forming stars, some disk material will enter the 
classical Hill volume and return again to the circumstellar 
disk. 

Second, at progressively lower background tempera- 
tures, the amplitudes of the spin variations increase, par- 
ticularly for the smaller enclosing sphere (bottom panels) 
where the variation increases by nearly a factor of two as 
measured by the rms deviations from the averages. The am- 
plitude and spatial extent of density variations in the enve- 
lope shows a similar trend, as we show in figure [TBI Within a 
distance of ~ l/4i?Hfrom the core, amplitudes of a factor of 
up 2-3 in normalized "peak to peak" variation and 30-70% 
in variance are present in all four simulations. At larger dis- 
tances, significant differences appear between the behaviors 
seen at different temperatures. While the variations in both 
amplitude and spatial extent decrease with increasing tem- 
perature in all four simulations, for the highest temperature 
realization they become essentially indistinguishable from 
the background flow within a distance of about 2_Rh ■ At the 
other extreme, for the lowest temperature realization, varia- 
tions remain significant up to ~ 4J?h from the core, and are 
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Figure 16. Variations in density as a function of distance from 
the core at t f» 74 yr (i.e. at the same time presented in figures [2} 
0, for simulations with different background temperatures, using 
the density extrema (top) and the density variance (bottom) as 
metrics, both normalized to the average density at each distance. 
The long dashed, solid, dotted, and short dashed lines define the 
variations for simulations tmJ^O, tm20, tmlO and tm05, respec- 
tively. The heavy solid line defines the variations intrinsic to the 
background flow, which are derived from the radial gradient of 
density as a function of distance from the star. 



non-zero all the way to the edge of the simulation volume, at 
6-Rh- The behavior of the middle temperature realizations 
are bracketed by those at the high and low extremes. 

We attribute the increases in variability, both in spin 
and in density, to feedback generated by perturbed mate- 
rial as it returns to the envelope for one or more additional 
encounters. Expelled material will inevitably exhibit differ- 
ences from the smooth background flow as it enters the core's 
environment. When temperatures are low, flow perturba- 
tions extend to larger distances and reach greater ampli- 
tudes than at higher temperatures, due to the decreased im- 
portance of pressure derived forces far from the core. These 
forces act to restore the flow to its background state and, 
when they are small, larger excursions may be sustained be- 
fore the restoring forces reach magnitudes large enough to 
return the flow to its unperturbed state. Significantly, vari- 
ability is not suppressed for any background temperature 
studied, and we conclude that the activity is robust. 



4-4-3 Generalizing the results to other core masses 

The simulations in the last section consider only IOMq 
cores. The results, however, can be generalized by framing 
them in terms of the the accretion radius, to parameterize 
the spatial extent of the activity, and the dimensionless ratio 
of the accretion radius to the Hill radius: 

Ra = G&M.) 1 ' 3 M*{ 3 

Rh dpi cf 

to parameterize characteristics of the flow activity. The ra- 
tio is inversely proportional to the background temperature, 
by way of the squared sound speed term, c 2 s , and its rela- 
tionship to temperature through the equation of state. It is 
also directly related to the core's mass present in the defi- 
nitions of both Ra and Rh, resulting in the proportionality 
of M p i 2//3 . Therefore, for a given semi-major axis, we can 
specify characteristics of the dynamical activity for any pair 
of (T, M p i) values. 

We can interpret the meaning of a given ratio by re- 
calling the physical meanings of the two quantities that de- 
fine it. Specifically, the accretion radius, as noted in section 
12.1.61 defines the location at which gas's thermal energy and 
gravitational energy of the core become comparable in mag- 
nitude. In turn, the Hill radius (or more generally, the full 
three dimensional Roche surface, for which the Hill radius 
defines the largest separation from the core) defines the loca- 
tions at which the centrifugal potential energy arising from 
the rotating frame and the combined gravitational potential 
energies of the star and core become comparable. Their ra- 
tio, therefore, becomes a measure of the relative importance 
of the rotating frame in comparison to thermodynamic prop- 
erties on the flow. 

In this context, we designate the region of parameter 
space where the ratio Ra/Ru is below unity, as the "ther- 
modynamic flow regime" , because the thermal energy of the 
background flow exceeds that of both gravitational and cen- 
trifugal potentials at the Hill radius. Therefore, the flow 
will be dominated by the pressure effects that define the 
thermodynamic state of the gas, with only small influences 
from the rotating frame. At smaller separations, effects of 
the rotating reference frame will have still smaller influence, 
while effects due to the gas's thermal properties increase. As 
gas falls deeper into the core's gravitational potential well, 
hydrodynamic feedback on the flow becomes important, as 
gravitational potential energy is converted into thermal en- 
ergy. We therefore expect flow feature s similar to the "non - 
stationary" behaviors reported by, e.g.. lRuffertl l|l997l . ll999l ). 
who report flow patterns around a point mass in which the 
direction of rotation oscillates in time. 

We designate the region of parameter space where the 
ratio is larger than unity, as the "celestial dynamic flow 
regime", because thermal energy plays a relatively smaller 
role at the Hill radius, so that gas flow more closely re- 
sembles that of a purely dynamical system of massless test 
particles. Systems in this regime will be characterized by 
prograde rotational flows near the core, a "horseshoe" or- 
bit region at distances outside the separatrix defined by the 
Hill volume, and circulating orbits both inside and outside of 
the orbit position of the core. Finally, for ratios near unity, 
thermal energies and gravitational energies, along with the 
forces associated with each, are of comparable magnitudes. 
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Figure 17. Lines of constant Ra/Ru, plotted in (T, M p i an ) 
space, together with the loci of several of our simulations, as 
marked. Three lines are plotted, showing ratios of unity (solid) 
and ratios where R\ is four times larger (dotted) or smaller 
(dashed) than ijjj. The regimes designated as celestial dynamic 
and thermodynamic flow regimes, as defined in the text, are also 
shown. 



We therefore may expect the flow to exhibit features of both 
types of flow. 

Figure [17] displays three lines of constant radius ratio 
in (T, Mpi) space, each assuming semi-major axis of 5.2 AU, 
as in our simulations. Following the discussion above, we ex- 
pect flow patterns of similar character at each point along 
these lines, and different flow patterns as the ratio changes. 
Specifically, that oscillating flows are ubiquitous for all of the 
simulations, as expected from their placement in the ther- 
modynamic flow regime. In addition, for lower temperatures, 
where the ratio lies closer to unity, the effects of the rotat- 
ing frame play a larger and larger role as the time averaged 
spin tends toward positive, or at least less negative, values. 
Higher amplitude oscillations are also generated, consistent 
with the increased importance of the celestial dynamic ef- 
fects due to the rotating frame. 

Of note in the figure is that for the low core masses typ- 
ical of the early stages of growth, and for temperatures of 
< 150 — 200 K, as expected in the region of core formation 
(i.e. the "ice line" , where solid densities increase due to the 
formation of ices), the flow characteristics fall well within 
the thermodynamic flow regime. Therefore, given a back- 
ground temperature near that of ice formation, we expect 
that a 1M® core will exhibit similar flow characteristics to 
our simulation tm^O, with an oscillating flow pattern present 
in a limited volume close to the core, with low amplitude. 
At the same temperature, a 3-5M® core would generate ac- 
tivity over a larger fraction of its Hill volume, with similar 
behavior to that of our prototype simulation tm20. Thermo- 
dynamic driving of the flow will become unimportant only 
at very high core masses or very low temperatures. 

We conclude that thermodynamically active flows will 
be an important characteristic of core environments during 
their early stages of growth. The activity will only decrease 
in importance for objects above a few tens of earth masses in 
size, as values of the radius ratio increase to well above unity. 



Because we have not attempted simulations with such high 
core masses however, we cannot determine with certainty 
the masses for which the activity will become negligible. 



5 RELEVANCE OF OUR MODELS FOR 
CHONDRULE FORMATION 

An important implication of the temperature and density 
data derived from our simulations, discussed in section 13.2.1 1 
lies well outside of the study of dynamical phenomena that 
may develop in the Jovian planet formation environment. 
Specifically, the densities, temperatures and time scales that 
occur in the envelope may be important for the theory of 
meteoritics and, more specifically, to models of chondrule 
formation. 

Chondrules, and the meteorites that contain them, rep- 
resent a large fraction of the volume of meteorites that 
have fallen to the earth and been studied. They are solidi- 
fied, spherical melt droplets of refractory materials typically 
up to a few milli meters in size. Briefly summarized from 
Ijones et al\ l|2000l ). we note for our purposes that chondrules 
underwent both a very rapid heating event to ~2000 K, 
followed quickly by a rapid cooling event of magnitude 50- 
1000 K/hr. Among a wide variety of models purporting to 
produce such conditions, passage of solid material through 
nebular shocks is currently favor ed as among the most likely 
(see e.g. lDesch fc Connolly|2002i ). Among this model's draw- 
backs are, first, that shocks that have the right density, tem- 
perature and velocity characteristics are hard to form and, 
second, that it is difficult to arrange for these shocks to exist 
for a long enough time to produce enough bodies to match 
the current observations. 

The parameter space for which chondrule production 
may occur in shocks is bounded by pre-shock densities 
within a factor of a few of 10 -9 g cm -3 , temperatures near 
300 K and shock speeds near 6-7 km s _1 (|Desch fc Connolly! 
2002, Table 5). For these conditions, the authors also note 
that enhancements in the particle density were important 
for fo r mation mo d els. C oncurrent work of iCiesla fc Hoodl 
(|2002T) : llida et al\ (|200ll ) come to similar conclusions. Fig- 
ures [2] and [4] show that appropriate background conditions 
exist in our simulations, and we propose that dynamical ac- 
tivity in the Jovian envelope could provide a source for both 
shocks and reversible compressive heating that remedy the 
shortcomings noted above. 

Although the basic temperature and density conditions 
that occur in our simulations can be seen, e.g., in figures [3] 
and[6j ascertaining whether short duration heating and cool- 
ing events are also present requires additional analysis. We 
considered and discarded the option of including test par- 
ticles that could be passively advected with the local flow, 
because of the significant added computational and storage 
cost they would require. Instead, we rely on the similar-but 
not identical-solution of 'advecting' test particles through a 
fixed-time snapshot of the conditions at a specific time. 

To that end, we have performed a streamline analysis 
of the trajectories of an ensemble of test particles injected 
into the simulation restart dump of our prototype model, 
for the same time as that shown for figures [5] and 2] Parti- 
cles are placed at a set of locations at the edge of the grid 
and allowed to advect through the simulation volume us- 
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Figure 18. The temperature (left panels) and density (right panels) encountered by three test particles as they are advected through a 
snapshot of the simulation volume, each as functions of a fictitious time coordinate, described in the text. The zero point for time has 
been arbitrarily shifted in each case so that the peak temperature occurs at ~ 10 days. 



ing the local fluid velocity, linearly interpolated between the 
grid zones adjacent to the particle at each time. Each par- 
ticle is given a time step based on that used to advance the 
gas at that location, so that it advances forwards through 
a fictitious time coordinate, passing through the volume at 
the rate of matter in the local flow. Similar linear interpo- 



lations of density and internal energy are used to derive the 
temperature from the equation of state. 

In Figure 1181 we show temperatures and densities for 
three test particles for which a passage through the envi- 
ronment of the core has occurred. The particles shown were 
chosen to illustrate the range of peak temperatures and den- 
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sities that may be encountered by slightly different trajec- 
tories through the envelope. Each encountered a short du- 
ration heating and cooling event as they passed through the 
dynamically active region close to the core, but the mag- 
nitudes and duration varied in each case. In the top ex- 
ample (test particle 1), the peak temperature rose to well 
over 3000 K and the density to nearly 10~ 8 g cm -3 as the 
particle's trajectory passed through the innermost regions 
of the envelope. The conditions encountered by the particle 
2 were much more moderate, with peak temperature and 
density values of ~ 1800 K and 4 x 10~ 9 g cm~ 3 respec- 
tively. Particle 3's trajectory took it only through the outer 
portion of the envelope, so that it encountered only much 
lower temperatures and densities, although in this case for 
a much longer period of time than either of the other two 
cases shown. All three particles encountered several days of 
cooler processing near 500-800 K, and a visual scan of many 
other similar events shows that such additional annealing is 
common. 

The temperature peaks for test particles 1 and 2 offer 
widths of < 1 day, with both a very rapid rise and fall. Close 
examination of their trajectories reveal that the widths re- 
flect essentially the crossing time for a single grid zone in the 
simulation. Therefore, although already quite narrow, they 
are actually likely to be overestimates of the true widths 
that would be obtained from the same models realized at 
still higher resolution. The dynamical activity in the en- 
velope, coupled with the temporally very narrow temper- 
ature/density peaks, especially in cases similar to that of 
particle 2, offer evidence that chondrule production could 
occur in the environment of Jovian planet formation. 



6 DISCUSSION 

6.1 Summary of our main results 

We have performed a suite of 3D numerical simulations of a 
cubic 'cutout' region around a 10M® point mass embedded 
in a circumstellar disk, in a modified shearing box frame- 
work. We find that the flow is not characterized by a simple 
one dimensional circumplanetary envelope embedded in a 
near static background flow. Instead it is highly dynamic, 
varying on time scales ranging from months to years far 
from the core, down to hours or days, very close to it. 

We used the RMS amplitude of time variations observed 
in components of specific angular momentum of material 
flowing around the core as a quantitative proxy for this ac- 
tivity and studied its sensitivity to changes in physical and 
numerical parameters used to model the evolution. Specifi- 
cally, we varied the character of the core, via its mass and 
gravitational softening coefficient, and we studied the sensi- 
tivity of the activity to artificial changes in the background 
flow, imposed by changes in the overall temperature or den- 
sity gradients in the disk, imposed by changes in the back- 
ground temperature of the disk, or imposed by changes in 
the assumed equation of state. We found that the ampli- 
tude is sensitive to the character of the core, increasing with 
core mass and decreasing with larger gravitational soften- 
ing. Increases in the background temperature of the disk 
cause decreases in both the amplitude and spatial extent of 
the activity. Activity can also nearly be suppressed when 



we assume a fixed equation of state-either isothermal or 
isentropic-which in turn serves as a proxy for more complex 
heating and cooling processes neglected in our current work. 
Each of these factors contribute to the efficiency and char- 
acter of hydrodynamic feedback on the large scale flow by 
the material returning to the background after encountering 
the core and we conclude that the origin of the activity is 
closely linked to such feedback. 

Finally, we demonstrate that an important consequence 
of the dynamical activity is to generate very short duration 
heating and cooling events in material that travels through 
the envelope. We compare conditions in these events to those 
expected to be required to form chondrules, and show that 
the ranges of temperatures and densities reached span those 
for which chondrule formation is expected. They also extend 
to longer duration and lower temperature events for which 
significant annealing may occur, and we propose that the 
origin Jovian planets and the origin of of chondrules and 
other annealed solids in the solar nebula are linked, through 
the dynamical activity in the envelope. If such a link can be 
established, then chondrules and annealed solids represent a 
physically examinable record of the processes present during 
the formation of Jovian planets. 

6.2 Potential implications of our results 

Our models represent a significant step towards a complete 
picture of the formation of Jovian planets in the context 
of the core accretion model. All of our simulations generate 
some dynamical activity and only simulations evolved using 
fixed EOS evolution or large gravitational softening produce 
activity of significantly lower magnitude than in our proto- 
type simulation and its many variations. Given the severe 
and unphysical restrictions on the thermodynamics imposed 
by fixed EOS evolution, we do not believe that the results 
of those simulations closely resemble the actual evolution of 
any real system. A large gravitational softening parameter is 
similarly restrictive, because it implies a spatially extended 
and relatively massive envelope, at very early phases of the 
core's growth. We therefore believe that more sophisticated 
models will show that dynamical activity will be present at 
levels important enough to play a role in overall formation 
process of Jovian planets. If so, our results mean that the 
standard core accretion model for Jovian planet formation 
will require revisions, to include the effects of such dynami- 
cal activity. 

Given this statement, we now speculate on what those 
effects may be, and their contribution to the planet forma- 
tion process. Most importantly, we expect the mass accre- 
tion rate to differ from that derived in the context of the 
standard core accretion model, because in that model the 
accretion is throttled by the rate of energy loss from the 
outer surface of the envelope. When the rate is low, pressure 
forces build up and choke off further accretion, when it is 
high, accretion proceeds more quickly. However, if the enve- 
lope has no stable outer surface, high pressures at any point 
on the 'surface' (such that one can be defined at all) will not 
hinder material flow through an adjacent region where they 
are low. Indeed, such high and low pressure regions may sim- 
ply be consequences of material flowing into or out of the 
core's environment. They may also be consequences of ther- 
mal energy generation from shocks created by interacting 
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material flows. In this case, an accurate model of the evolu- 
tion will require an accurate accounting of the contributions 
from such energy generation and advective terms, to both 
the mass accretion and energy loss rates of the material in 
the core's environment. 

Depending on its strength, dynamical activity may ei- 
ther delay accretion, by repeatedly disrupting the formation 
of the envelope, or enhance it, by driving instabilities in the 
flow. In the context of our current models, we cannot deter- 
mine which outcome is more likely, or if both occur at dif- 
ferent times. If the envelope is disrupted on timescales short 
compared to its cooling timescale, then our present mod- 
els, which neglect cooling entirely, will provide an accurate 
picture of the long term flow pattern. Any material which 
does cool will be swept away, rather than becoming a perma- 
nent part of the core's envelope. If instead instabilities occur 
which do not entirely disrupt the envelope, a likely conse- 
quence will be rapid mass and angular momentum trans- 
port through the envelope onto the co re. Instabilities are 
known to develop in rotating polytropes (|Pickett et al .1 19961 ; 
llbnian et al\ \l998). morphologically similar to the planet 
formation environments discussed here, but whether those 
or other, similar instabilities can be triggered in the present 
context awaits further study. Also, the consequences of such 
rapid accretion are unclear. How quickly would a massive 
envelope develop in which further activity did not play a 
large role? 

The scenario we present offers a number of attractive 
advantages over other models of Jovian planet formation. 
Unlike both the class of formation models using global grav- 
itational instabilities in the circumstellar disk and of many 
previous models of core accretion, a massive disk is not re- 
quired at any time during the evolution. In particular, no 
massive disk is required to exist for the long period of time 
required for the classic picture of core accretion, for which 
most gas accretion occurs only at the conclusion of a long 
hydrostatic growth period. Assuming that the activity re- 
sults in enhanced mass accretion rates, it also provides a 
formation mechanism which avoids one of the most critical 
unsolved issues in standard models of core accretion, that 
the core does not grow quickly enough to avoid migrating 
through the disk and onto the star. Finally, our scenario 
naturally provides a mechanism for producing very short 
duration heating and cooling events with thermodynamic 
characteristics similar to those expected to produce chon- 
drules and annealed solids. Production in this scenario will 
endure for a significant fraction of the formation time scale 
of Jovian planets (itself a significant fraction of the disk life- 
time), resulting both in a large yield of objects and allowing 
both processing and reprocessing events to occur, both fac- 
tors consistent with the meteorite record. 



6.3 Comparisons to other work 

The work described in this paper can be profitably 
compared to three main categories of research on Jo- 
vian planets. Together they make up a spectrum of 
models ranging on one extreme, to multidimensional 
models which explore the interaction of a circumstellar 
disk with a massive object embedd ed within it, using 
relatively simple physical m odels dLubow et all 1 19991 ; 
iD'Angelo etaU 120021 . l2003bl lal; iBate et all 120031 ). and on 



the other extreme, to one dimensional models which 
explore the core accre tion paradigm using highly de- 
tailed physical models 
1 2004 iHubickvi efal 
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attention is paid to the interactions between the envelope 
and the surrounding disk. 

The main focus of most previous multi-dimensional 
work has been on large scale interactions between the core 
and the disk as a whole, in order to study the migration of 
the core through the disk and determine parameters related 
to its survival. In this context, many details of the struc- 
ture of the envelope become less important, and approxi- 
mations are made to sidestep the more difficult numerical 
and physical issues. In comparison, our focus has been pri- 
marily on the interaction of the envelope with the core and 
the immediately surrounding disk material. Accordingly, we 
provide only a rudimentary description of the background 
disk, through a boundary condition and through its gravi- 
tational effect on the material inside our simulation volume. 
We concentrate our spatial resolution on the envelope re- 
gion only, thereby permitting much finer resolution of the 
flow features. The work presented here represents the high- 
est resolution simulations studying the formation of Jovian 
planets so far published. Our work also includes a physical 
model which, while necessarily more simple than required 
for a complete picture of the evolution, at least avoids some 
of the most restrictive assumptions of previous models. 

As noted in section |3~T1 the morphologies found in our 
work differ from those of the previous studies, exhibiting 
both more dynamical activity and less disk-like or spiral 
structures in the circum-core environment. The reasons for 
the differences between the morphology most likely lie in dif- 
ferences in the physical assumptions made in each. Specifi- 
cally, the most important assumptions include those regard- 
ing the thermodynamic treatment for the gas, those regard- 
ing its distribution around the core and those regarding the 
continued mass accretion onto the core itself. We now dis- 
cuss the physical implications of each of these assumptions. 

We employ an ideal gas equation of state to close the 
system of hydrodynamic equations. Additionally, no sources 
of heating or cooling are included apart from shock heat- 
ing generated by the dynamical interactions of the flow on 
itself. In contrast, many previous workers have chosen an 
isothermal equation of state, which artificially (and very ef- 
ficiently) removes thermal energy generated by those same 
interactions. High temperatures and pressures are therefore 
suppressed, and the feedback loop that would otherwise ex- 
ist as the hot gas perturbs incoming material cannot form. 
As we showed in section [4.4.1l our simulations do not exhibit 
dynamical activity when an isothermal equation of state is 
used, consistent with the lack of such activity in other work. 
A complete model will of course include some cooling effects 
due to radiation, which will undoubtedly act to suppress 
the activity we see, however the magnitude of this effect is 
difficult to estimate from the currently available literature. 

Although some calculations have been performed 
(AB09) which include radiative cooling, they do not super- 
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sede our results due to limitations implied by the resolution 
afforded by the simulations themselves. Their spatial reso- 
lution, although high in relative terms, is still much lower 
than is required resolve the flow sufficiently at small dis- 
tances from the core to capture the extremely small scale 
features that characterize the shocks and dynamical activ- 
ity close to the core. For example, tests of the shock cap- 
turing characteristics of the PPM algorithm in use in our 
simulatons, show that shocks can be resolved with a maxi- 
mum of 2-3 grid zones, corresponding to distances of < 2Rj 
at the grid zoning used in our simulations. In contrast, the 
SPH ('Smoothed Particle Hydrodynamics') method used by 
AB09 requires at least 5-10 particl e smoothing lengths to re- 
solve a shock (see, e.g., section 6 of Price 2012). While AB09 
quote a minimum smoothing length over their simulations 
comparable to our zone size, their overall particle count is 
quite low, implying that very few particles have such small 
sizes, while most are much larger. Coupled with the com- 
paratively large artificial viscosity required by the method, 
we believe that the features we observe in our simulations 
would not have been observable in theirs, independent of 
any considerations regarding cooling. 

Our work assumes that no gas has accumulated around 
the core. This assumption manifests explicitly in our ini- 
tial conditions, for which we define only a 'bare' core em- 
bedded in an unperturbed background flow, and implicitly 
through our specification of the gravitational softening of 
the core. In addition to its normal purpose as an aid to 
sidestep various numerical issues in finitely resolved simula- 
tions, using non-zero gravitational softening is equivalent to 
the assumption that the core's mass is distributed over a fi- 
nite spatial volum e, corresponding to the s patial scale of the 
softening (see, e.g.. |P'Angelo et al\ (|2003al ) for a particularly 
detailed discussion of such softened potentials). The gravi- 
tational softening used our simulations is smaller than em- 
ployed in many other works, where simulations with value s 
of ~i? H /10-i? H /5 are typical (e.g. iD'Angelo et all l2003al ). 
In contrast, the largest softening we use in any simulation 
is (for simulation so 10), while in all others at the 

same resolution it is a factor 16 smaller. Given the results 
of section 14.21 where we showed that softening values this 
large are sufficient to suppress dynamical activity, we con- 
clude that other published work would not have seen such 
activity, due to the differences between their assumptions 
and ours. 

Further, the effect of permitting the core to accrete es- 
sentially all material moving into a small volume around 
it is to suppress any buildup of gas in the envelope. Such 
buildup would otherwise provide a back pressure and either 
stall further accretion or, as happens here, provides a sup- 
ply of material with which additional infalling material can 
interact. For simulations of cores with sufficient mass, well 
abov e that expected to c ause the onset of 'core instabil- 
ity' ijWuchterl e t al. 2000), unimpeded accretion may be a 
physically relevant and important process to consider. In the 
present case however, where we study the evolution around 
only 10M® cores with very low mass envelopes, such an as- 
sumption is not valid. Unimpeded accretion is also physically 
inconsistent with the existence of the massive envelope im- 
plied by the large gravitational softening parameters, since 
any such envelope would quickly be accreted into the core 
itself. 



Our spatial resolution is nearly fine enough to resolve 
the core (see section [2.1.5l) . In this context, implementing a 
'black hole' accretion model, in which all material moving 
into a predefined volume around the core is accreted, would 
be inconsistent with the known character of the core (i.e. 
that it is solid). Such material becomes unavailable to inter- 
act with other material falling onto the core at later times. 
Therefore, models in which accretion is permitted may unin- 
tentially suppress the feedback loop that generates the vig- 
orous dynamical activity we observe in our simulations. 

Previous ID work focuses primarily on the secular 
growth of the envelope and core over its long term (10 5 -10 8 
year) evolution. It therefore includes detailed treatments of 
physical processes relevant over long time scales, while as- 
suming that short term dynamical fluctuations average out 
over time. In comparison, our models include a simpler phys- 
ical model but fully 3D spatial resolution. Due to the compu- 
tational cost of 3D models however, our simulations extend 
over an very short period of the total evolutionary history 
of the planet's growth. 

Each of these ID models still requires that a bound- 
ary condition be applied at the outer edge of the enve- 
lope. To the meager extent that defining an envelope can 
be meaningful in the simulations we present, our dynam- 
ical results ar e broadly cons i stent with the condition as- 
sumed in the iLissauer et all l|20Qg| ) simulations, in which 
the boun dary was defined at a radius of ~7?h/4 from the 
core. The lLissauer et al\ (|2009h assumption is based on the 
result of dynamical calculat ions somewhat similar to our 
own (|D'Angelo et al\ 1200 3a - ), in which only tracer parti- 
cles closer than 7?h/4 from the core remained near it. Our 
simulations do not form any envelope structure, and mate- 
rial in the volume that they define as most strongly bound 
to the core is most dynamically active in our simulations. 
We do however, concur with the result that material at 
larger separations (still within the Hill volume) is largely 
composed of background disk material whose trajectory has 
merely caused it to pass near the core, rather than become 
bound to it. This same resul t is inconsistent with that of 
iTerquem fc Heinemannl l|201ll ) , who find that the volume of 
gas that is bound to the core will grow to fill its entire Roche 
lobe even if it initially fills only some fraction of it. 

An important concern regarding all ID models is the 
issue of whether or not the assumption that hydrodynamic 
properties of the system actually do average out over time, 
so that a hydrostatic evolutionary model is accurate. Do hy- 
drodynamic effects force the evolution to proceed in a direc- 
tion that it otherwise would not have taken, considering only 
the hydrostatic properties of the system? Our 3D models il- 
lustrate that the formation environment for Jovian planets 
cannot be adequately described by any hydrostatic treat- 
ment. Such a treatment would imply that material and en- 
ergy transfer between the envelope and disk occurs through 
a stable boundary, with material far from that boundary 
being only indirectly affected. On the other hand, our simu- 
lations show that the interactions between disk and envelope 
frequently disrupt the envelope entirely, so that its growth 
must begin again from scratch. Also, with such vigorous 
activity, cooling by any process becomes far less relevant, 
since the cooled material does not remain near the core for 
long enough to affect the long term evolution of the enve- 
lope. Since our simulations extend over only a tiny fraction 
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of the total formation history of Jovian planets and also 
do not incorporate a sufficiently complete physical model of 
the evolution, we cannot yet speculate as to the fraction of 
a planet's formation history such a statement will be true. 
In the following two sections, we discuss the implications of 
such activity on current core accretion models, and the com- 
ponents of physical models that we believe will be required 
to make such conclusions. 



7 UNANSWERED QUESTIONS 

The physical model included in our work is a necessarily 
simplistic view of the actual environment of forming Jovian 
planets and, because of that, our work generates more ques- 
tions than it answers about the growth of Jovian planets 
and solar system formation. In this section, we point out a 
number of such simplifications in our models, discuss what 
we believe is the relative importance of each in the context 
of Jovian planet formation, and pose a number of questions 
that can be addressed as they are removed. 



7.1 Simplifications that affect the thermodynamic 
state of the envelope 

We made two simplifications in this work which appear to 
be critical to making definitive connections between our re- 
sults and the physical systems they are intended to model. 
Each of these factors contribute directly to models of the 
internal structure of a hydrostatic (or nearly so) envelope, 
determining its energy balance and its density and temper- 
ature structure as a function of radius. They are critical 
components of both hydrodynamic and hydrostatic models 
in this context. Before any detailed analysis of the conditions 
will be of lasting value for theories of planet formation, each 
must be accurately accounted for, particularly in the deeper 
regions of the envelope where dynamical activity is most 
vigorous. 

First, no account has been taken of the radiative cooling 
and heating processes that undoubtedly occur in the enve- 
lope, beyond our models employing fixed equations of state. 
As discussed in the context of the p ossibilities for grav i- 
tational instabilities to form planets (|Durisen et al\ [20071 ). 
use of fixed equations of state as proxies for cooling models 
greatly restricts the behavior of the energy balance in the 
gas. Whether or not these restrictions correctly represent the 
actual heating and cooling behavior of the Jovian envelope 
and its environment, or instead artificially suppress the ac- 
tivity which should actually be present, must be addressed 
in follow-up work. 

In our simulations, the envelope volume is character- 
ized both by high temperatures and densities and by rapid 
changes and large gradients of the same quantities. It seems 
likely that including the effects of radiative transfer will tend 
to smooth out temperature gradients throughout the sys- 
tem. How much so, and to what extent are the temperatures 
and densities artificially altered from their actual values by 
our physical model? To what extent will the dynamical prop- 
erties of the system will change when they are included? Will 
the dynamical activity continue? Will shocks still develop in 
the flow? Preliminary indications with the comparatively 



extreme restrictions of locally isothermal and locally isen- 
tropic equations of state suggest that at least some activity 
remains, so we remain hopeful that our current treatment is 
approximately correct. 

Unfortunately, the material opacities needed to model 
radiative transport accurately are poorly constrained, since 
they arise largely from the properties of ice and dust, both in 
terms of grain sizes distribution and in terms of their compo- 
sition. In many regions throughout the envelope, solids may 
be vaporized either temporarily or permanently, causing 
opacities to change by orders of magnitude over small spa- 
tial and temporal scales. When solids are removed as opacity 
sources, the envelope may become optically thin so that ra- 
diative cooling becomes more efficient, particularly closer to 
the core where temperatures are highest. Will radiative cool- 
ing in these regions redistribute or remove enough energy 
from the gas to change the activity that occurs when radia- 
tive processes are neglected? Moreover, when they reform, 
grains with similar chemical compositions may have quite 
different opacities than originally. Radiative transfer mod- 
els which incorporate this grain history will likely require 
the advection of tracers in the simulations, which model dif- 
ferent chemical compositions of different solid species. This 
is challenging, both in the context of developing a physi- 
cal model and numerically. To what extent will such models 
affect the dynamics? 

Second, no account has been taken of the temperature 
and density dependence of the equation of state. Including a 
complete equation of state will also cause important changes 
to the flow because, at various temperatures, the adiabatic 
exponent, 7, of the gas will differ from the constant value of 
7 = 1.42 that we have assumed. While this value will be ap- 
proximately correct for temperatures relevant for the back- 
ground circumstellar disk material (i.e. T ~ 100 — 1000 K), 
where at various temperatures rotational and vibrational 
modes of hydrogen molecules play a greater or lesser role, 
it becomes much less so for higher temperatures and densi- 
ties characteristic of the gas deep in the envelope near the 
core. Above T ~ 1500 — 2000 K, dissociation and ionization 
become increasingly important, and the effective adiabatic 
exponent falls to values as low as 7 ~ 1.1, before increasing 
again to 7 = 5/3 when the gas becomes fully dissociated. 
As shown in figure [4] the envelope temperatures span this 
range. For a given compression, a lower adiabatic exponent 
means that pressure increases by a smaller amount. There- 
fore, our fixed exponent treatment may artificially increase 
pressure gradients beyond the physically relevant levels they 
should reach, thereby artificially stimulating dynamical ac- 
tivity to a greater than realistic extent. On the other hand, 
lower exponents mean that a given compression event can 
increase the mass density to a higher values, perhaps to lev- 
els at which self gravitating instabilities can develop. Which 
overall effect will dominate? 

Once these thermodynamic effects are included, it will 
be of interest to run simulations for far longer in time, in 
order to explore questions relating to whether activity con- 
tinues indefinitely or instead decays with time. Our current 
models simulate only a 100 yr segment of the Jovian planet 
formation history, during a time when the envelope did not 
contain much mass. Because there were no energy loss mech- 
anisms in our current models, no mass accretion into the en- 
velope could be expected. Will dynamical activity continue 
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to overwhelm net mass accretion onto the core, will activ- 
ity continue while superposed upon an overall net accretion 
rate (and at what rate?), or will the activity simply decay 
over time to levels unimportant for the evolution? 

Over such long timescales, A third simplification of 
some significance to the thermodynamics will also play a 
role. Namely, the effect of heating due to planetesimals in- 
falling into the envelope. Over the short periods simulated in 
our models, such heating is not likely to play a large role in 
the dynamical state, because only a few events occur dur- 
ing such a short time and because the energy input from 
any single event is not a large fraction of the internal en- 
ergy of the envelope. Such conditions will not be true over 
longer periods, during which many more events occur. How 
will dynamical activity change when this additional heating 
source is accounted for? In our current simulations, disk ma- 
terial enters the envelope then leaves it again, carrying with 
it thermal energy generated in shocks or other interactions 
with material near the core. Will a similar outcome hold for 
the heating generated by an infalling planetesimal, or will 
the additional thermal energy be trapped in the envelope, 
escaping only slowly to the background disk? 

Longer simulations will permit studies not only of the 
accretion rate, but also of the envelope's behavior as it be- 
comes more massive. Does activity depend on having only a 
low mass envelope? Does activity enhance the net accretion 
rate or reduce it? An answer to this question is important in 
the context of the overall survival of a proto-planet because 
the observed lifetimes of the circumstellar disks out of which 
plane t form may only be 2-4 mill ion year s (Haiscf Tet all 
2001), while core accretion models (|Pollack et aZ.lll996T ) re- 
quire as many as ~6 million years even in a relatively mas- 
sive massive disks (2-3 times the minimum solar nebula). 

7.2 Simplifications in the initial conditions and 
the physical realization of the disk and core 

In addition to the simplifications in the thermodynamic 
treatment, our models also study only a limited range of 
the parameter space of interest in terms of the initial condi- 
tions explored. They employ simplified treatments of several 
physically relevant components of the overall system, such 
as the treatment of the core itself, and the vertical structure 
of the disk. 

In this work, we have studied the flow around only 
IOMq cores, but such objects must first, of course, grow from 
smaller masses. Will dynamical activity occur in the flow 
around lower mass cores? The discussion in section 14.4.31 
would suggests that activity will occur, but would be less 
vigorous for the smallest cores for conditions typical in the 
parts of the disk where cores are expected to form. Is there 
a lower limit below which activity does not play an impor- 
tant role in the evolution? If the character of mass accretion 
also changes markedly with the strength of the dynamical 
activity, then an important observable consequence of such a 
limit would be the present day masses of the Jovian planets 
in our solar system. Can such a correlation be made? 

Due largely to limits imposed by computational cost, 
we have approximated the core as a softened point mass 
with no actual size or surface. With only a factor of ~ 2 in- 
crease in spatial resolution from the highest employed here, 
simulations would begin to resolve the core and some, more 



physical, accounting must then be made of its structure. 
We forsee that a simple solid surface realized by a reflecting 
boundary condition will be the next step, and will provide 
a physical picture adequate to model the influence of the 
core's spatial extent on the flow. Of interest will be the ques- 
tion of whether or not interactions with the surface enhance 
dynamical activity, perhaps by being reflected off of it to in- 
teract with other infalling material. Also, the gravitational 
potential well reaches its maximum depth at the core's sur- 
face. Resolving the flow there means that the temperatures 
and densities seen close to the core will represent the most 
extreme of those expected in any dynamical flow. Will some 
fraction of such material also be ejected into the background 
disk flow? Will it carry signatures of its passage through the 
core environment? 

A still more complex treatment of the surface would 
require models of the interaction between the material in 
the core and the flow. To what extent does mass exchange 
between the core and envelope occur? Does high-Z core ma- 
terial mix into the envelope, (and perhaps also back out 
into the disk itself?), thereby enriching it, or does it remain 
largely unaffected by the envelope activity? Does such mix- 
ing change with envelope mass? 

Indirectly, the accretion rate will also influence the 
migration rate of the core through the disk. As its mass 
increases, the mutual gravitational torques between the 
core and disk will increase in turn and the core inter- 
acts more strongly. Whether the stronger interaction leads 
to increased migration will depend on details of the mass 
distribution around the core. Accou nting only for torques 
from Lindblad resonances (see e.g., I Ward! 1 19973 ), migra- 
tion will accelerate. Accounting for corotation torques as 
well, migration may decelerate or even reverse direction 
l|Peplinski. Artvmowicz fc Mellemair2008al fbT). 

While omitting gravitational forces in the z coordinate 
precludes a description of the disk's vertical structure, we 
expect that the influence of that structure on the envelope 
activity will not be large. Here again, the activity is limited 
only to the immediate environment of the core, a factor of 
several in spatial extent smaller than the disk's scale height. 
A possibly more important consequence of such a treatment 
will be the possibility of mass being ejected entirely out of 
the disk, to fall later on some region distant to the core. In 
this case, the question of how large the accretor must be 
before it first begins to deplete mass from an entire vertical 
column of the disk. For an accretor of that size, accretion 
may be reduced due to the loss of material accreting onto 
the poles of the proto-planet and a gap may form due to 
the accretion of most of the material near the planet, also 
slowing additional migration through the disk. 

Of somewhat lesser importance for the near term are the 
geometric approximations imposed by our choice of coordi- 
nate system and our omission of background gravitational 
forces in the z (vertical) coordinate of our disks. Mitigat- 
ing the importance of both cases, is the fact that activity is 
limited to a comparatively small volume close to the core. 
Therefore, the coordinate system itself plays only a small 
role. 
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7.3 Is there a connection between models of 

Jovian planet formation models and models of 
chondrule formation? 

As noted above, a consequence of dynamical activity in the 
envelope is that material is not irrevocably bound to the core 
after first entering the envelope environment. Instead, some 
fraction of the material returns to the circumstellar disk 
from which it originated. Will this material retain a signa- 
ture of its passage through the hot, circumplanetary region? 
The possibility is intriguing because if so, some such mate- 
rial could still exist today, potentially providing a record of 
the conditions present during the formation of Jovian plan- 
ets that could be directly observed. In fa ct, a common class 
of m eteorites-chondrites-exist (see e.g. iPalme fc Bovntonl 
1 19931 ) and are found throughout the solar system, which 
exhibit signatures of high temperature processing on time 
scales of only a few hours or days. Can the formation of 
Jupiter actually be linked to the formation of chondrules ob- 
served in the meteorite record though a model such as ours? 
Can it be linked to the formation of othe r materials, such 
as th e annealed silicates found in comets (|Harker &: Deschl 
2002), for which the required temperatures and densities are 
lower? 

Although the results from our initial streamline analy- 
sis are promising, they are no substitute for an investigation 
of the trajectories of specific packets of material through the 
system. In order to establish any real connection, we must 
perform a much more detailed analysis of the conditions en- 
countered by such packets. Is a particle's thermodynamic 
trajectory the same when it is advected through an actual 
time dependent flow, as opposed to the fictitious advection 
through a fixed flow that we have performed? Will more de- 
tailed analysis show that the shocks that are generated fit 
into the required density/temperature/velocity parameter 
space? One concern already apparent is that the flow ve- 
locities of material flowing through the shocks (1-2 km s _1 , 
as estimated from the directly available fluid flow veloci- 
ties thems elves) are uncomfo r tably low compared to thos e 
quoted by iDesch fc Connolly! (|2002fi and llida et ail (|200lf ). 
It seems unlikely that the velocities will be increased as dra- 
matically as that by any of the improvements to the models 
we might make. Finally, what fraction of material encoun- 
ters conditions appropriate for chondrule formation during 
its passage, in comparison to material that instead encoun- 
ters regions that are inappropriate? And what fraction of the 
total budget of solid material in the solar nebula undergoes 
such processing? Equivalently, do conditions appropriate for 
chondrule formation exist for a large or small fraction of the 
disk lifetime, so that a larger or smaller fraction of solids 
encounter those conditions during their evolution? Finally, 
how do processed materials get from where they form (near 
5 AU) to their final locations, in meteorites throughout the 
inner solar system? 



7.4 Future Directions 

We believe the most profitable course following this work 
will be to investigate the character of activity that is gener- 
ated for cores of different masses, in order to confirm directly 
the scaling relationship discussed in section P4.4.3I In paral- 
lel with this effort, models employing an accurate equation 



of state for the circumstellar material will be of great in- 
terest. Such an equation of state must include the various 
molecular and atomic states of hydrogen, which will likely 
cause changes to the strength of the hydrodynamic feed- 
back through the changes they imply for the ratio of specific 
heats, 7. Ultimately however, we expect that models includ- 
ing radiative transfer will be required to accurately model 
the full process of mass flow and accretion in the neighbor- 
hood of forming planetary cores. In future work, we hope to 
explore all of these questions. 



APPENDIX A: DIFFERENCES FROM THE 
STANDARD SHEARING SHEET FORMALISM 

The differences between the standard shearing sheet model 
and our modifications to it are to be found in the equa- 
tions of motion in the x and y directions, equations [3] and 
[4] corresponding to an underlying radial and azimuth coor- 
dinate frame that is rotating with the angular speed of the 
core. Additional differences are found in the specification of 
the boundary treatments, as discussed in section [2~3l above. 
To illustrate the differences in mathematical formalisms we 
show here the approximations underlying both our formal- 
ism and the shearing sheet and their differences, starting 
from the equations of motion in a non-rotating cylindrical 
coordinate frame. 

In an inertial, cylindrical coordinate frame, the equation 
of motion in the radial and azimuthal directions are 

d(pVr) | d{pV r V r ) | d{pV T V„,) | d(pV r V z ) pV$ 

dt dr rd(f> dz r 

= -|-P^ (AD 
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where the capital letters on the velocities are used to denote 
the inertial frame velocities. In eauation lAll the pV^/r term 
defines an analogue of the centrifugal force which, mathe- 
matically, comes from the fact that the full time derivative 
of the momentum includes components that account for the 
fact that direction vectors f and <f> are dependent on position 
and through them, also time. In equation IA2I the pV r V$/r 
term plays a similar role, as a near analogue of the Coriolis 
force. 

When translating to a frame rotating at angular velocity 
Slfr, we must include centrifugal and Coriolis terms of the 
form fif r x (fifr x r) — 2fif r x v. Additionally, the meanings 
of each of the variables as defined above changes, to reflect 
the moving coordinate system. For purposes of clarity, we 
note the identifications (V r — ¥ V r , V$ — > rfl{ T + v$, V z — V v z ) 
for each of the coordinate direction. In other words, lower 
case variable names denote quantities in the moving frame. 
Then, the form of equations IA1I and IA2I change to include 
the additional terms as follows: 
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and 

d(pv<f,) djpv^vr) djj>v$v$) djpVtfyVz) pv r v^ 
dt dr rd<f> dz r 

dv 

+2pv r Q b = £■ - p—. (A4) 

ro(f> raq> 

The cross term in equation IA3I arises from the Coriolis 
term, 2f2/ r x v, while the first and third terms come from 
the centrifugal forces due to the curvi-linear coordinate sys- 
tem and the rotating frame, respectively. Together, the three 
terms in eqution IA3I form a perfect square identical to the 
corresponding single term in equation I All We complete the 
specification of the equations of motion defined in equations 
[3] and [4] by replacing the cylindrical coordinate variables 
(r, r<fi) with (a;, y). 

In both the treatment above and the shearing sheet, 
the identification of the coordinates pairs (x,r) and (y,r<f>) 
is made in each of the hydrodynamic equations. The shear- 
ing sheet approximation goes further than this by omitting 
the Coriolis-like term due to the curvature of the coordinate 
system (i.e. the term pv r v^/r, above) and expands the sum 
of the centrifugal force and the gravitational force of the 
central star in a Taylor series, for which the zeroth order 
term is assumed to be zero, and only the linear term is re- 
tained. lHawlev. Gammie fc Balbusl (l995h state the form of 
the equations to be solved in the shearing sheet approxima- 
tion, cast as partial derivatives of velocities. For comparison 
to the form used here, we recast these equations in momen- 
tum form, omitting all terms relating to magnetic fields, 
which we neglect entirely: 

djpvx) d(pv x v x ) d(pv x v y ) d{pv x v z ) 
dt dx dy dz 

-2pv y Q {r = + 2pqV?x (A5) 
ox 

O(pVy) , d(pVyV X ) d(pVyVy) d{pV y V Z ) 

dt dx dy dz 

+2pv x Q b =~ (A6) 
dy 

where q — — d In Q^/dln r ~ 3/2, for a nearly Keplerian 
disk. The last term on the right hand side of equation I A5I 
defines the first order approximation to the effective poten- 
tial in which the gas flows noted above. Rather than use 
this Taylor expansion approach, we retain the full form of 
the background potential of the star and disk, assuming only 
that it is constant in time. This permits a direct specifica- 
tion of the gravitational potential to be made, rather than 
an analytic approximation. It also permits transient veloc- 
ity fluctuations that temporarily upset the assumption that 
centrifugal and gravitational forces are in equilibrium to be 
accounted for more fully. 



APPENDIX B: NUMERICAL ISSUES 
REGARDING TREATMENT OF THE 
BOUNDARIES 

The simplest method of producing a flow-out boundary in 
either the inner or outer x ('radial') direction would be to 



reproduce the flow variables of the last interior grid zone 
in the boundary cells adjacent to that zone and just out- 
side the grid. We have found this method unsatisfactory for 
our simulations because the background state is effectively a 
dynamical equilibrium condition rather than a static equilib- 
rium condition. Implementing the simplest flow-out bound- 
ary changes the gradients there (e.g. the pressure gradient 
partially responsible for the underlying rotation curve), and 
is equivalent to violating the dynamical equilibrium condi- 
tion. The effect of an initially small outflow at the negative 
x boundary, corresponding to the inner radial boundary of 
the underlying disk, is to flatten the gradient and further 
strengthen the outflow. At the positive x boundary a flat- 
tened gradient has the opposite effect: the effect of an ini- 
tially small outflow is to temporarily create a very strong, 
local inflow condition. 

Since the major gradients of the background flow are in 
the radial direction (x direction in the simulation cube), the 
same concern does not apply to the y boundaries, however 
they can still affected by another problem. If the disparity 
between the smooth background state and the flow near the 
boundary is large, the simulated flow can still be adversely 
affected. In particular, the fixed condition at the boundary 
means that if some volume of gas flows towards that bound- 
ary but does not approximately match the conditions there, 
it may either be accelerated as it leaves the grid, causing 
mass further inside the simulation volume to be similarly 
accelerated into the prematurely evacuated volume, or it 
may be decelerated, effectively 'piling up' at the boundary 
and perhaps reentering the flow. 

We have closely monitored the flow in our simulations 
for consequences of this effect and have observed it in only 
one (labeled tm05 in table[T]above). We believe that the fail- 
ure of this simulation is due primarily to the fact that the 
assumed background temperature was very low. As a conse- 
quence the Jeans wavelength of the mass in the grid became 
roughly comparable to the linear dimension of the simula- 
tion cube itself. As the planet drew additional mass into its 
influence, the perturbation on the local potential grew as 
well and further amplified the effect. The conditions just in- 
side the boundary were then no longer similar to those just 
outside it and violated our smooth background assumption, 
resulting in matter turning back into the simulation volume 
and perturbing the evolution still further. A separate sim- 
ulation, neglecting the local disk self gravity, did not suffer 
from this defect, and we discuss it in sect ion |4 . 4 . 2 1 b elow as 
a more physically realistic alternative to simulation tm05, 
with a nearly identical physical model to its higher temper- 
ature cousins. 
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